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ABSTRACT 


The  time  histories  of  emission  lines  from  succes- 
sive ionization  stages  of  impurity  atoms  in  a plasma 
are  predicted  by  solving  the  coupled  rate  equations 
governing  the  populations  of  these  ionization  stages 
with  the  help  of  a computer  program.  A detailed  de- 
scription of  the  program  is  given.  An  application  of 
such  predictions  for  the  measurement  of  ionization 
rate  coefficients  of  ions  occuring  in  a rapidly 
ionizing  plasma  is  discussed. 


I . INTRODUCTION 


Knowledge  of  the  lifetime  of  various  ionic  stages  of  impurity 
atoms  in  high  temperature  hydrogen  or  deuterium  plasmas  is  very  im- 
portant in  controlled  nuclear  fusion  research  for  the  estimation  of 
energy  losses  as  well  as  in  understanding,  e.g.,  the  physics  of  the 
solar  corona  (flares,  etc.).  The  basic  processes  governing  the 
ionic  lifetimes  are  collisional  ionization  and  the  various  recom- 
bination processes — radiative,  dielectronic,  three-bodv,  whose  rela- 
tive importance  is  determined  by  the  plasma  conditions.  In  order 
to  understand  the  relative  importance  of  these  processes,  the  plasma 
conditions  could  be  simulated  and  the  time  histories  of  emission 
lines  from  various  ionization  stages  be  predicted  by  a computer  pro- 
gram. The  purpose  of  this  report  is  to  describe  in  detail  such  a 

computer  code.^  It  was  used  for  example,  to  determine  the  ioniza- 
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tion  rate  coefficients  of  various  ions  using  a ©-pinch  plasma  * 

with  added  impurities.  The  principles  of  such  measurements  are  as 

a 

described  by  Hinnov. 

In  Section  II  the  basic  equations  and  the  adopted  procedure 
for  their  solution  are  described.  In  Section  III  the  set  up  of  the 
program  is  described.  In  Section  IV  the  description  and  listing 
of  the  program  is  given.  In  Section  V sample  test  data  are  set 
up  such  that  analytical  solutions  can  be  obtained.  A comparison 
with  the  computer  solution  checks  the  accuracy.  In  Appendix  I 

plotting  subroutines  are  listed.  In  Appendix  II  an  application  of 
this  program  for  the  determination  of  ionization  rates  is  discussed 
using  the  data  from  an  experiment^  as  an  example. 


II.  BASIC  EQUATIONS 


The  impurity  ions  in  a transient  plasma  go  successively  through 
the  various  ionization  stages,  the  degree  of  ionization  lagging  be- 
hind the  corresponding  quasistationary  corona  equilibrium. 

The  density  N^  of  ions  (K  is  the  nuclear  charge)  in  the  total  plasma 
Volume  V is  determined  by  the  following  rate  equation: 


d <NkV) 
dt 


VN,  . NI,  , - VN,  NI.  + 
k-1  k-1  k k 


VN.  , Na,  , - VN,  Na. 
k+1  k+1  k k 


(1) 


where  N is  the  electron  density,  1^  is  the  ionization  rate  coeffi- 
cient and  a^  the  recombination  rate  coefficient.  The  plasma  volume 
would  change  with  time  in  a 0-pinch  plasma  due  to  compression  or  ex- 
pansion. This  variation  was  accounted  for  in  our  experiments  by  the  as- 
sumption that  particle  losses  are  negligible  , which  is  found  to  be 
valid  in  the  operation  of  our  0-pinch  with  a reverse  bias  field. 

d (NV)  - 0. 

dt 

dV  V dN 

dt  " “ N dt  (2) 

This  also  assumes  small  concentrations  of  impurity  ions  such  that  elec- 
tron production  fcy  ionization  is  negligible.  Hence  using  Eq.  2 in  1 


Sc 

dt 


- N. 


k-1 


NIk-r  \ NIk + 


Vi  N<Vi  - \Nak + ± \ 

N at 


(3) 

However,  for  different  plasma  machines  and  plasma  conditions 
Eq.  (2)  could  be  different  which  would  make  the  last  term  in  Eq.  (3) 


s 


to  be  different.  Hence  the  last  term  in  Eq.  (3)  Is  called  the  source 
term  which  should  be  appropriately  modeled. 

These  coupled  rate  equations  are  solved  numerically  with  the 
measured  electron  density  and  temperature  given  as  Input. 

The  following  formulae  were  used  for  the  rate  coefficients  In 
Eq.  (3): 

Ionization:  A choice  can  be  made  to  use  one  of  the  following 

approximations . 

2 

1.  Semiempirical  formula  due  to  Kunze  derived  by  fiting  to 
the  semlemperlcal  predictions  of  Lotz.8  It  is  referred  to 
here  as  "S(k)-Lotz  & Kunze" 


S(K) 


SA(K) 


(*■>  x) 3 + 40 


kT*  -E./kT 
Ei+3kT  6 


3 

cm 


/sec 


(A) 


SA(K)  - 7.5  x 10~ 8 ni 

Ei  (5) 

Here  "th."  is  the  number  of  electrons  in  the  1 subshell  (e.g. 

2 2 2 
for  the  Is  2s  2p  configuration  is  1 for  2p,  2 for  2s 

2 

and  2 for  Is  subshells,  and  the  ionization  rate  of  the 

ground  state  is  a sum  of  all  these  contributions)  and  E^ 

is  the  corresponding  ionization  energy.  Both  kT  and  E^ 

are  in  e.v.  The  constant  SA(K)  is  given  as  input. 

A comparison  to  another  semlemperlcal  formula  due  to 
4 

Hinnov  was  found  to  yield  the  same  results  as  Eq.  (4). 

2.  Summers  described  the  procedure  for  calculating  the 

ionization  rate  coefficients  using  the  Exchange  Collisional 

- 3 - 


8 

and  Impact  Parameter  (E.C.I.P.)  method  due  to  Burgess. 

This  procedure  is  adopted  into  our  program  and  the  ioniza- 
tion rate  is  referred  to  as  S(K)-  Burgess  and  Summers. 

3 

S(K)  = A(K)  X (Burgess  and  Summers  rate.)  cm  /sec 

(6) 

A(K)  is  a multiplicative  value  to  change  the  rate  if  one 
wishes.  It  is  given  as  input.  Other  input  data  needed  to 
evaluate  Eq.  (6)  Is  described  in  the  Comments  of  Subroutine 
Main  - R/W  in  Sec.  IV. 

Recombination:  Three  processes  of  recombination  are  considered. 
Generally,  for  a rapidly  ionizing  plasma  recombination  rates  are 
negligible.  However,  one  has  to  check  this  especially  for  higher 
ionization  stages.  The  total  recombination  rate  is  referred  to 
as  Alpha (K) . 

1.  Radiative  recombination:  RALPHA 

For  the  evaluation  of  RALPHA  a 

9 

formula  due  to  Seaton  (Eq.  4.3.35  Ref.  10)  was 
adopted. 

RALPHA  (Rct^)  = 5.20.10"14  Z(E/kT)1/2  (0.429  + 0.5  Hn  + 

0.469  (kT/E)*^  (cm^sec  ^)  (7) 

where  'Z  ' is  the  effective  charge  of  the  recombining  ion  K (giving  K-l)  and 
E is  the  ionization  potential  of  the  ion  upon  recombination,  (i.e.  K-l) 
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Dielectronic  recombination:  DALPHA 


For  the  evaluation  of  DALPHA  simple  formulae  due  to  Lan- 
dini  and  Monaignori  Fossi11  are  used. 

For  ions  pertaining  to  H,  He,  Ne,  K to  Ni  DALPHA  (Da  ) = 

k 

1.60  x 10~10  (Z+l)2  F(K)(kT)~3/2  WOQ'^xp  [~---9^34W^K)j 

(8) 

where  Z is  the  charge  of  the  recombining  ion,  W(K) 
is  the  first  excitation  energy  in  e.v  for  the  recombining 
ion  and  F(K)  could  be  taken  as  the  number  of  electrons  in  the 
outer  shell  (F(k)=i  for  H-iike,  2 for  He-like,  6 for  Ne-like 
...) 

For  ions  pertaining  to  Li  to  F,  NatoAr,  CutoKr  DALPHA  (D&  ) = 

k 

1.6  x 10~10  (Z+l)2  F(K)(kT)*exp  £ ~0*9134  ^ 

(8) 

where  Z+l  and  W(K)  are  defined  as  in  Eq.  8«  F (K)  values 
are  taken  from  Ref.  7. 

Three  body  recombination:  BALPHA 

For  the  evaluation  of  three  body  recombination  a formula 
„ , 12 

due  to  Grlem  was  used. 

BALPHA  (b^)  - 1.4  x 10*31  N Z~6 (e±  /klj  *e xp  [%  / kTCn^l)2^ 

(9) 


♦ 


where  N is  the  electron  density,  Z the  charge  of  the  ion 

before  recombination  and  n1  the  quantum  number  of  the  Colli- 
12 

sion  limit  i.e.  the  level  from  which  radiative  decay  is 
about  as  probable  as  collisional  excitation  to  higher  levels. 


i 2 14/17  -2/17  -1/17 

s 1.26  x 1<T  Z N (E./kT)  ^ exp 


4 EiA7(n1)3  kT 


(10) 


A plot  of  n1  is  given  in  Ref.  9. 

Having  determined  the  time  history  of  ionic  populations, 
the  time  histories  of  their  emission  lines  could  be  predicted. 
The  emission  coefficient  of  an  optically  thin  allowed  line, 
whose  upper  level  is  mostly  populated  by  electron  collisions 
from  the  ground  state  is  given  by 


€ = N X \ ergs  / cm3  sr  sec  (11) 

We  assume  excited  state  populations  to  be  small  so  that 

represents  the  ground-state  population*  Also  X is  the  excita- 
tion rate  coefficient.  If  N and  X are  constant  > the  time 

history  of  the  line  is  identical  to  the  time  history 
of  the  ion.  However,  variations  in  N and  X can  be  accounted 
for  even  though  the  absolute  value  of  X is  not  known.  The 
time  dependence  of  X is  accounted  for  by  the  effective  Gaunt- 
f actor  approximation. 


X » 1.6  x 10'5  f e"AE/kT  (AE) 1 (kT)'1/2 
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(12) 
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where  AE  is  the  excitation  potential  in  ev  and  ' f'  is  the 
oscillator  strength.  This  should  be  true  at  least  for  the 
usually  small  temperature  variations  occuring  during  the 
emission  of  many  lines  in  our  0-pinch  plasma.  However, 
such  an  assumption  has  to  be  supported  by  the  experimental 
determination  of  density  and  temperature  profiles. 

III.  SET  UP  OF  THE  PROGRAM 

The  code  has  a main  programs  four  subroutines  and  three 
external  functions.  They  are  described  below. 

1.  MAIN-R/W:  The  main  program  used  to  read  in  and  print  out 
the  data,  and  control  the  rest  of  the  program. 

2.  SOLVE:  A subroutine  that  solves  the  differential  equations. 

3.  RATES:  A subroutine  that  interpolates  density  and  tempera- 
ture, calculates  the  rate  coefficients  and  increments  in 
populations  at  each  fine  time  step  'J'  required  by  SOLVE. 

External  functions  SUMl  and  SUM2  are  for  the  evaluation  of  S(K)  - 
Burgess  and  Summers. 

4.  PRINT:  Print  is  used  to  output  the  calculations  of  ionic 
populations  calculated  in  "SOLVE",  Print  together  with 
external  function  EIOFX  calculate  the  intensities  of  spec- 
tral lines. 

5.  PLOT:  The  intensity  as  a function  of  time  is  stored  for 
teletype  plotting. 


The  program  TV  PLOT  and  the  Subroutines  TT  PLT  1,  TT  PLT  2,  and 
TT  PLT  3 are  for  teletype  plotting  the  data  stored  by  PLOT. 


IV.  PROGRAM  DETAILS  AND  LISTING 

The  following  variables  are  in  the  common  block  of  each  Sub- 
routine. 

J.  Index  counter  counting  the  fine  time  steps  in  the  Subroutine 
SOLVE.  When  J reaches  Jmax  the  Subroutine  PRINT  is  called  and 
the  output  is  printed.  Then  J is  set  back  to  3. 
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Jmax:  Given  as  data  in  Main-R/W.  The  maximum  value  J takes  is 

Jmax. 

ND:  Number  of  ions  ION)  for  which  SOLVE  solves  the  differen- 

tial  equations.  It  is  given  as  input. 

H:  Variable  representing  the  step  size (time)  in  SOLVE. 

DTmax:  Predicted  value  of  H given  as  input. 

Kttep:  Counter  counting  the  number  of  fine  steps  in  the  integra- 

tion program  SOLVE. 

KSend:  Maximum  limit  "Kstep"  could  take  in  the  Subroutine  SOLVE. 

It  is  given  as  input. 

KDOB:  Flag.  If  it  becomes  -1  the  time  step  H is  considered  for 

doubling  in  SOLVE  subroutine. 

KHAV1  Flags  used  in  SOLVE  to  decide  whether  at  a certain  time, the 
KHAV2 

KHAV3  step  size  H is  to  be  made  half. 

D0B1,  Counters  used  in  SOLVE  to  limit  the  maximum  number  of 
DOB  2 

iterations  in  finding  proper  H. 

T(6Q):  Time  at  each  fine  time  step  J in  SOLVE. 

EPSI(60):  Accumulated  error  indicating  how  well  the  integration  is 
proceeding. 

Y(60,18):  Population  density  of  each  ionization  stage  at  each  time 
step. 

YR(60,18):  Increments  in  populations  for  each  time  step  — i.e. 

These  are  calculated  in  "RATES". 

S(18):  Ionization  rate  Coefficients. 

Alpha(18) : Total  recombination  rate  Coefficients. 


NTIME: 


Maximum  number  of  coarse  time  steps  for  the  input  data. 


TIME  (100),  DENS  (100),  TEMP  (100):  Electron  density  and  temperature 
data  at  coarse  time  steps  L=1  to  NTIME  given  as  input. 

DNE(60),  TE(60) : Interpolated  electron  density  and  temperature  data 
at  each  five  time  step  J. 

Ion:  Maximum  number  of  ions  for  which  rate  equations  are  going  to  be 
solved. 

IonMl:  Ion  -1. 

KI0N(18) : Counter  for  the  number  of  ions. 

SA(18) : The  coefficients  for  the  ionization  rates  calculated  by 

the  formula  of  Lotz  and  Kunze.  They  are  given  as  in- 
put. 

SB(18) : Values  used  in  Lotz  and  Kunze  formula.  (SB  is  set  3.) 

SE(18):  Ionization  potentials  for  the  ions  given  as  input. 

ANPRIM(18) : Data  for  three  body  recombination  in  RATES  to  be  given 

as  input. 

NAME(18):  Identification  of  the  plasma  condition. 

SFN(18):  Oscillator  strengths  for  the  lines  whose  time  histories 

are  to  be  predicted. 

SEX(18) : Excitation  potentials  in  eV  for  these  lines. 

001:  Electron  density  divided  by  the  sum  of  all  ionic  popu- 

lations at  t=0  ysecs. 

-4 

YMIN:  10  times  the  sum  of  all  ionic  populations. 

YYMIN:  10  ^ times  the  sum  of  all  ionic  populations. 

EN,  ED,  EH:  Constants  to  check  accuracy  in  SOLVE. 

NDT,  NAINI:  Constants  used  in  the  program. 

NA:  A variable  representing  the  stage  of  ionization  below 

which  populations  are  less  than  YMIN. 
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Input 


THI(60,  18): 


F(18) : 
W(18) : 

Yes (18) : 

Flag: 

Tim: 

DIF  2: 

TTT(60) : 

Flag  1: 


M(I,K): 

H(I,K): 


A(K)  : 


Theoretical  intensities  stored  for  plotting, 
data  for  dielectronic  recombination  (DALPHA)  formulae 
due  to  Landini  and  Fossi.^ 

Given  in  Ref.  11. 

First  excitation  energy  in  recombining  ion.  Given  in 
Ref.  11. 

A flag  whose  value  is  read  in  for  each  ion  to  decide 
which  formula  for  DALPHA  is  to  be  used. 

A flag  decides  whether  experimental  or  theoretical  in- 
tensities are  to  be  plotted. 

Variable  used  in  "PRINT"  to  find  the  fine  time  steps 
that  correspond  to  the  Coarse  time  steps. 

The  difference  between  fine  time  step  and  the  coarse 
time  step  evaluated  in  "PRINT". 

The  fine  time  steps  that  are  equal  to  or  close  to  the 
coarse  time  steps. 

A flag  read  in  input  data  which  decides  whether  Burgess 
and  Summers  or  Lotz  and  Kunze  ionization  rate  coeffi- 
cient is  to  be  used. 

Number  of  equivalent  electrons  in  shell  I of  ion  K. 
Ionization  energy  in  eV  of  an  electron  in  shell  "I"  of 
ion  "K".  Five  shells  should  be  more  than  enough. 

A factor  multiplying  Burgess  and  Summers  ionization 
rate  to  match  time  histories. 


PROGRAM  DETAILS: 

1.  MAIN  - R/W:  It  reads  in  the  data  and  manages  the  rest  of  sub- 
routines. The  integer  variables  FLAG,  FLAG1  and  YES  are  flags  whose 
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values  given  as  input  would  make  the  required  calculations.  The  pos- 
sible manipulations  are  described  in  the  comments  of  the  listing  of 
this  subroutine.  In  the  appendix  II  data  is  set  up  as  an  example. 

The  array  YI(J,K)  is  used  to  read  in  the  experimental  time  histories. 

A brief  description  of  the  purpose  of  the  following  variables  is  given 
below. 

A SUM:  Sum  of  initial  populations  of  all  ions. 

001:  Electron  density  at  time  t=0  us  divided  by  A SUM. 

This  is  used  in  PRINT  to  normalize  the  calculated 
populations.  This  procedure  is  to  check  how  the 
effect  of  compression  was  taken  care  of  in  the  sub- 
routine RATES  by  the  variable  ABL  and  how  the 
smooth  transfer  of  populations  from  one  stage  to 
the  other  is  proceeding. 

YMin  and  YYMin:  These  are  used  to  avoid  the  accuracy  check  in  SOLVE 

for  ions  with  negligible  populations  at  a certain 
time  step  J. 

2.  SOLVE:  This  subroutine  solves  the  coupled  differential  equations 
by  numerical  integration  using  Adams  predictor-corrector  methods. 

A general  outline  of  the  working  process  of  this  subroutine  is  de- 
scribed below. 

The  flag  NSTART  takes  values  1,  2 or  3 given  as  input  data.  If 
NSTART  is  set  to  1,  the  program,  between  statement  numbers  5 through 
19,  will  try  to  find  a proper  time  step  H for  continuing  the  solution 
with  stability  and  convergence.  An  initial  value  for  H (DTMAX) 

is  given  as  input.  Also  the  density  and  temperature  information  has 


to  be  given  at  one  negative  time  (non  physical)  for  this  self  starting 
procedure.  This  information  is  only  useful  for  starting  the  program 
and  does  not  influence  the  solution. 

If  NSTART  is  set  to  2,  the  initial  value  for  H given  as  input  is 
used  to  continue  the  integration.  However,  the  program  tries  to  adjust 
"H"  to  obtain  solutions  within  the  accuracy  requirements. 

If  NSTART  is  set  to  3,  the  difference  between  coarse  time  steps 
1 and  2 at  which  density  and  temperature  information  was  given  would 
be  taken  as  the  starting  value  of  H and  the  program  tries  to  adjust 

H,  if  required,  to  continue  integration  keeping  the  required  accuracy. 

The  program  had  provisions  to  speed  up  the  calculation  by  doubling 
the  step  size  or  to  retrace  back  one  time  step  and  cut  the  step  size 
to  half  and  continue  the  integration  within  the  accuracy  requirements. 

The  accuracy  checks  were  done  by  using  the  values  for  EN,  EH  and  ED  as 

I,  .land  .001  respectively.  These  values  were  found  to  yield  good 
stability,  convergence  and  accuracy  in  the  running  of  SOLVE.  In  sec- 
tion V analytical  solutions  were  compared  with  numerical  solutions 
for  a simple  two  ion  system.  Further  details  about  the  working  of  the 
program  are  in  the  comments  of  the  listing  of  SOLVE. 

3.  RATES:  It  calculates  the  ionization  and  recombination  rates,  as 

well  as  increments  in  populations  at  each  time  step  J.  The  variable  ABL  takes 
care  of  effects  of  any  compression  etc.  The  external  functions  SUM1  and 
SUM2  are  used  to  calculate  S(K) — Burgess  and  Summers  if  desired 
(i.e.  Flag  1.GT.0).  Other  details  are  in  comments  of  the  listing. 
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PRINT:  This  subroutine  prints  the  concentration  of  each  ionization 

stage  at  each  time  step  J.  It  calculates  the  intensities  of  emmission 
lines,  using  the  external  function  EIOFX  and  prints  them. 


MAIN  - R/V 


1.  C 

2.  C 

3.  C 

4.  C 

5.  C 

6.  C 

7.  C 

8.  C 

9.  C 
10. 

11. 

12. 

13. 

14. 

15. 

16. 

17. 

18. 

19. 

20. 

21.  C 

22.  C 

2*3 . c 

24.  C 

23.  C 

26.  C 

27.  C 

28.  C 

29.  C 

30.  C 

31.  C 

32.  C 

33.  C 

34.  C 

35.  C 


THIS  PROGRAM  IS  FOR  SOLVING  COUPLED  RATE  EQUATIONS  AND  PREDICTING 
TIME  HISTORIES  OF  THE  IMPURITT  IONS  AND  THEIR  EMISSON  LINES  IN  A 
RAPIDLY  10N12ING  PLASMA  LIKE  THAT  OF  A THETA  PINCH. 

DATA  NEEDED  AS  INPUT  IS  ELECTRON  DENSITY  AND  TEMPERATURE  PROFILES 
WITH  RESPECT  TO  TIME  AND  THE  C H A RE C TE R I S TI C DATA  FOR  EVALUATING 
THE  VARIOUS  RATE  C 0 E F F I C I EN T S--I ON Z AT I ON &R E C OMB IN AT  10 N (R AD I AT  I V E , 
THREE  BODYR 01  ELECTRONIC ) RATE  COEFFICIENTS. 


INTE6ER  TIM, FLAG, FLAGl , YES 
DIMENSION  YI (60,18) 

COMMON  J ,JMAX,ND ,H,DTMAX , ND T , K S T EP , KS E ND , E ND , E D , E H , EN , KO OB , 

1 K FAV1 ,KHAV2 ,KHAV3 , DOB  1 ,D0B2 

COMMON  T (60, EPSIC60),  Y ( 60 , 1 8 ) , YR  ( 60 , 1 8 
1 >,ANPRIM(1S) 

COMMON  S (1 8) , ALPHA (18 ) ,DENS (100) ,TEMP<100) .TIME  <1 00) ,SA( 18) ,SB (18) 
1 , SE  (18) .NAME (8) ,K I ON (18) , DNE (60) ,TE ( 60) , SF N ( 1 8 ) , S E X ( 18 ) 

COMMON  NT IME , I ON , I0NM1 ,NA,NAINI ,001 ,YMIN,YYMIN 
COMMON  TH(6  0,18),TIM,FLAG,F(16),YES(18),W(18),M(5,18),HI(5,18) 

1 ,FLAG1, TTT(60),DIF2,A(18) 


•« 

.# 


IN  THIS  PROGRAM  COARSE  TIME  STEPS  MEAN  THE  TIME  STEPS  AT  NHICH 
DENSITY  AND  TEMPERATURE  INFORMATION  IS  GIVEN. 

TIM  -IT  IS  A VARIABLE  USED  IN  THE  SUBROUTINE  PRINT  TO  FIND  THE 

COARSE  TIME  STEPS  AMONG  THE  FINE  TIME  STEPS  GENERATED  BY  THE 
SUBROUTINE  SOLVE. 

FLAG  IS  A FLAG  TO  INSTRUCT  THE  COMPUTER  WHAT  IS  TO  BE  PLOTTED. 

'FLAG'  VALLES: 

0 - EXPERIMENTAL  INTENSITIES  ONLY  (PLOT  THEM) 

1 - BOTF  THEORETICAL  AND  EXPERIMENTAL  INTENSITIES  (PLOT  BOTH 

OF  THEM.  (THE  THEORETICAL  INTENSITIES  ARE  THE  ONES  COMPUTED  BY 
THE  SUBROUTINE  SOLVE.)  THIS  VALUE  CAN  ONLY  BE  USED  FOR 
REMOTE  TERMINAL  PLOTTING  BUT  NOT  FOR  TELETYPE  PLOT 


■■i 

j 

■j 

5 


I 


\ 

I 


36.  c this  is  the  way  sub  routines  ttstore  s plotter  are  construc 

3?.  C -TED. 

38.  C 2 - THEORETICAL  INTENSITIES  ONLY  (PLOT  THEM) 

39.  C 

40.  C FLA61  IS  A FLAG  TO  INSTUCT  THE  COMPUTER  WHICH  FORMULAS  FOR  IONIZATION 

41.  C ARE  TO  BE  U SED. 

42.  C FLAG1  VALUES: 

43.  C 0 - SEMI  EMPERICAL  FORMULA  DUE  TO  LOTZRKUNZE  IS  USED. 

44.  C 1 - SEMI  CLASSICAL  FORMULAE  DUE  TO  BU RG E S S & S UMM E RS  ARE 

45.  C USED. TWO  SUBROUTINES  S UM 1 & SUM  2 DO  THIS  JOB. 

46.  C EXPLANATION  OF  THE  INPUT  DATA  FOP  THE  RATE  COEFFICIENTS  IS  GIVEN 

47.  C BELOW. 

48.  C SACK)  SSB(K)  ARE  CONSTANTS  USED  IN  LOTZKKUNZE'S  IONIZATION  RATE 

49.  C COEFF  IC IEN  T.  IF  FLAG"!  IS  NOT  SET  C THEN  THESE  ARE  STILL  READ  BUT 

5C.  C NEVER  USED  . 

51.  C 

52.  C THE  FOLLOWING  IS  FOR  SUMME R S &b U R GE S S IONIZATION  RATE  COEFFICIENT. 

53.  C ,,M(I»K)"-  NO  OF  EQUIVALENT  ELECTRONS  IN  "IMSHELL  OF  ION  K. 

54.  C "H(I,K)"*  IONIZATION  ENERGY  IN  EV  OF  AN  ELECTRON  IN  SHELL  "I"  OF  ION 

55.  C K • FIVE  SHELLS  SHOULD  BE  MORE  THAN  ENOUGH  ! 

56.  C 


C DALPHA  IS  DIELECTRONIC  RECOMBINATION.  THE  FOLLOWING  ARE  NEEDED  TO 
C IT'S  EVALUATION. 

C 'F'  IS  INPUT  DATA  FOR  DALPHA. 

C 'W'  -FIRST  EX  ENERGY  OF  RECOMBINING  ION-  IN  DALPHA. 

C 

READ<5,1  CO)  (NAME(N) ,K  = 1 ,8) 

100  FORMAT (8 A4  ) 

1 READ(5,1C1)  ION, NT1ME.NAINI .FLAG, FLAG1 

101  FORMAT  <1  216) 

READ(5,1C2)(KI0N(K),SA(K),SB(K),SE(K),Y(2,K),ANPRIM(K),K=1,10N) 

102  FORMATU  12  ,4E1  2 .4  , F 1 2 . 4) 

READ(5,1 C3)  JMAX,ND,NDT,KSEND,END,ED,EH,  EN.DTMAX 

103  FORMAT(3  14  ,ie,5Efe.1> 

READ (5 , 1 C4 ) (TIME(K),TEKP<K),DENS(K),K*1,NTIME) 

104  FORMATC3E12.4) 

REAOC5.1C5)  (SFN(K),SEX(K),Ks1,ION) 

105  FORMAT ( 2 FI C. 5 ) 


\ r A D ( S # 1 C A HVESOO.FCk  ) W ( n ) , K = 1 , 1 0 \ ) 

! ' . i "St-AM  1 3,  f'  .r  ,fn.!0 

i f (FLAt  1.  •.  r . -)  iO  T n ’ / 

^Z'  k-1,I0\ 

•V  r A P ( 5 , 1 * " ) A ( K > , ( ( 1 , O , H I < I , Y ) , I = 1 , l ) 

1 ' F f:  K v a T ( H ' . ( 5 l ! , f ' .<■)>> 

if  (flac  ,(T.  i)  oorr  t 

\ T = \ T I ••  < - * 

D f'  11"  J = 1 , NT 
r TT  ( J ) = 7 IV:CJ+W) 

1 r A r ( r'  , 1 ’ '■  ) ( V I ( J , * ) , F - 1 , I . \ ) 

* F T < 1 :F  A.  2 ) 

/■Ll  'Hr  ‘ ■ F F r c L DATA  If.  READ  !Y  !.  C\,  . I T WAS  DOME  IN  THE  FOLLOWING 

''•'ft'.  DATA  F CP  < U N Z E N L 0 T C , P t N S J TV  i T E N P P R 0 F I l £ f,  , 0 A T A FOK  EXCITATION 

, ;•  AL’  H a , si!*  *F  1.  c i,Au9Et  ss  , « t X FEUKMAL  TIvt  HiST0PIES. 

C r 'L  r ' J = 1 , * T 

' «riTL(t.  ,11  Li  mU,K),>sI,10M 
1c  F r •• .<!  A T ( 1 c F f . 2 ) 

*.  * N ) F N A L I 2 r A ,\  p f*  L C T T H t EXt-fcPlNtNTAL  INTENSITIES 

CALL  P L 0 T T ( Y I , 1 , c ) 

IF  ( F L » L . l fv  . „ ) STOP 


- r 

i C A 

1 C 1 
< *< 


! 7 ’ ;ONt  1 = 10 r.  -1 

*P!TE  << , zc;> 

"C  f C»'-  A T < 1 FT  ,1  ( X #7  f.M  J 0 N I ? A ” I T N (' f INPU«1TY  IONS  IN  A PLASMA  WITH  61 
U r >4  DENSITY  AND  TENT'S®  A Tim  ///) 

:ci  FCRf'AT  (1  :x  ,1  :h  IN  PURITY  , A «.  / / ) 


I'?. 

«PJ  TE ( 6,  1) 

K-. 

• PITS ( 6 » Z't) 

1 C*  . 

- •: ; facpaki  cx,: 

i ir. 

if  (Flag  i • i n 

in. 

IE  (*,  7K) 

17  - 


I 

I 


i 


i 

i 


i 

I 


112. 
113. 
1 14. 
1 15. 
1 16. 
1 17. 
1 18. 
1 19. 
1 2C  • 
1 21  . 
1 22. 
1 23. 
1 24. 
1 25. 
1 26  . 
1 27. 
1 28. 
1 29. 
1 30. 
1 31. 
1 32. 
1 33. 
1 34. 
135. 
1 36. 
137. 
1 38. 
1 39. 
140. 
141  . 

142. 

143. 


144. 

145. 

146. 

147. 

148. 

k.  i 1*’. 


312  F ORMAT ( 1 CX ,4. l H DATA  FOR  ECIP  IONIZATION  RATE  COEFFICIENT  //) 

313  K =1  ,1  ON 

313  toRITE (6  , 31O(K10N<K),<M(I,K.)t  HI  < I , K ) , 1 = 1,5),  Y(2,K),ANPR1M(K),A(K)) 
3U  F0RMAT(5X,I2,5(I2,E10.4),E1C.4,F3.1,F10.6) 

•RITE (6, 21 5) 

315  FORMAT ( // /I  OX  , ' THE  IONIZATION  RATE  COEFFICIENTS  ARE  CALCULATED 

1BY  THE  F CRM)LA'//1uX,'S  = A * S E M I C L A S S I C A L RATE  COEFFICIENT  ' 
2///10X,  ' THE  FORMULA  USED  FOR  TOTAL  RECOMBINATION : ' / / 

3 1 0X , 8SH  ALPhA  = 5.2E-14*Z*SQRT(E (I-1)/T)  * (3.429  + 0.5*LN(E(I- 

4 1 ) / T ) ♦ C. 469* CBRT ( T/E ( I -1) ) + / 1 5 X , 7 5 H 1 . 4 E -3 1 * N * (NPRIME/Z)* 
5*6  * (E<  1-1) /T )**2  * EXP(E(I-1 )/ (T* (NPRIME*1)**?) ) »2H  ♦ , / 1 5 X , 

6'  DIELECTRONIC  R E C OM B I N A T 1 0 N ( L A N D IN  1 & F 0 S S I SOLAR  PHYSICS  2G.322P 
7 (1971 ) ) *7  // 

£11X,'FOR  THE  EXCITATION  WE  USE  G(EFF)  APPROXIMATION''  / / 1 0 X , * X=  1.5 
9SF-C5*FN*P*(EXP-EX/T)/(EX*SGRT(T) )'////) 

GO  TO  21 t 
211  WRITE(6,Z02) 

233  FORMAT (1  CX  , ' CHARACTERISTIC  DATA  FOR  THE  IONS  '//15X, 

1 4 H ION,5X,2H  A,10X,2H  B , 1 C X 7 H E (EV)  , 10X, 'INITIAL  D EN  S I T Y ' , 1 5 X , 

2 6HNPRIME//) 

WRITE  (6*  204)  <M0N(L)»SA(L),SB(L),SE(L),Y(2,L),ANPRIM(L),L  = 1,I0N) 
2C4  fOR*"AT(1!X,I3,1X,3t12.4,12XtE12.4,10XtFl2.4) 

WRI TE (6, 205) 

205  FORMAT ( ///10X,'  THE  IONIZATION  RATE  COEFFICIENTS  ARE  CALCULATED 
1 B Y THE  FCRFULA'//1GX,  ' S = A * ( ( LN ( 4 0* T / E ) * * 3 ) ♦ 4 0 ) * SQ R T ( T ) * E XP ( - E 
2/T)/(E*BT)'///10X,  ' the  formula  USED  FOR  TOTAL  recombination:'// 

3 1CX,  85H  ALPHA  = 5.2E-14*Z*SCiRT(E  (I-1)/T)  * (C.420  * 0.5*LN(E(I- 

411/T)  ♦ C.469*CBRT (T/ECI-1) ) ♦ / 1 5 X , 75 H 1 . 4 E -3 1 * N * (NPRIME/Z)* 

5*6  * (EC  1-1) / T > * * 2 * EXP (E ( 1-1 ) / (T* (NPR1ME + 1 ) **2) ) , 2H  *,/15X  , 

6'  DIELECTRONIC  R E C 0MB  I N A T ION ( L AN D IN  I 8 F 0 S S I -S 0 L A R P H Y . 2 0 ( 1 9 7 1 ) ) ' / / / 
711X,'F0R  THE  EXCITATION  WE  USE  G(EFF)  A P P R 0 X I M A T I ON ' / / 1 0 X , ' X = 1.5 
86E-05*FN*P*(EXP-EX/T)/(EX*SQRT(T))'////) 

3 1*6  WRITE(6,  208)  J M A X , N D « NDT  , KS  E ND  , E ND  , ED  , EH  , E N , D TM  AX 
208  FORMATd  CX  .'CONSTANTS  FOR  THE  INTEGRATION  PROGRAM  '// 

1 1 0 X , 7 H Jw  A X =I4,5X,5X,5H  ND  =I4,5X,6H  N D T =14, 

2 5 X , 5H  KSEND  =It//1UX,6H  END  =E8.1,5X,5H  ED  =E8.1,5X,5H  EH  =E8.1, 

3 5 X , 5 H EN  = E 8 • 1 ,5X,°H  DTMAX  =E8.4 //) 

WRITE  (6,  214)  . 


150. 

151. 

152. 
1 53. 
1 54. 
155. 
1 56. 
1 5?. 
1 52  . 
1 59. 
160. 
161  . 
1 62. 
1 63. 
1 64  . 
165. 
1 66. 
1 67  . 
168. 
1 69. 
1 7C  • 
1 71. 
172. 
1 73. 
1 74. 
1 75. 
1 76. 
1 77. 
1 78. 
1 79. 
1 80. 
181. 
1 82. 
183. 
1 84. 
185. 
1 86. 
187. 


214  FORM AT (/  // 10X , 36H  DATA  FOP  DIELECTRONIC  RECOMBINATION//, 

* 1 CX , 3H 10 N , 3X , 7HF0RMUL A , 3X  ,7H  F.0R.N,gX,6H  » ( E V ) / / ) 

-RITE (6,  21 5)  (KION(K)  »YES (K)  , F ( N ) ,W(K)  , k = 1 , ICN) 

215  FORMAT (1 5X , I 3 , 7X ,I3,7X , F3  .0 , 2 X , F 1 D . 5 ) 
wRI TE  <6 , 20d) 

206  fOP»AT(1  H jlCX.'nPUT  DENSITY  i TEMPERATURE  PROFILES  ARE:'/// 

1 1 5 X , 5 F TIME,7X,11H  TEMP  (EV)  ,16H  D EN S 1 T Y ( C M * * - 3 ) / / ) 

WRITE  <6,  20  7)  <TIME(K),TEMP(K),DENS(K)fk=1,NTI"E) 

2C7  FORMAT ( 1CX.3E12.4) 

WRITE (6,  21 2) 

212  FORMATC/  /I CX, ' DATA  FOR  THE  EXCITATION  OF  THE  IONS  '//15X,4H  ION 
1 , 5X,3H  F N, £X ,6HE  < ( EV ) / / ) 

WRITE  <6,  213)  (KION(K) , S F N ( K ) , S L X ( K ) ,K=1,ION) 

213  FORMAT(15X,I3,1XfF12.5,1X,Fl2.3) 

WRITE  (6,  210 

210  FORMAT ( 1 H ) 

C 

C 

ASUM  = C.O 
DO  30  1 k = 1 , ION 
3C1  ASUM  = ASUM  ♦ Y ( 2 , K ) 

001  = DENS(2)/A  SUM 
YMIN  = 1.0E~04*ASUM 
YYMIN  = 1.0E-06*ASUM 
T ( 2 ) = C.O 
NA  = MINI 
TIM  = 3 

CALL  SOLVE ( 1 ,NEND ) 

CALL  PRINT 
WRITE  <6  ,209)  NEND 

209  FORMAT ( 1 FO ,50X ,7H  NEND  =l5//10x,  'H U R R A Y »••»!* S UC C E S S F UL  RUN  !!!' 
IF (NEND  ,GT .1)  60  T0  50C 
C 

C THE  FOLLOWING  PART  IS  FOR  PLOTTING: 

C 

DO  370  K *1  ,ION 
DO  360  I *3  ,NTIME 


to 


1. 

18?.  3 7 0 CONTINUE 

189.  DO  271  I =3  , N T I M E 

19C.  371  TTTCI-2) =TTT(I) 

i 

191.  NTTrNTIME-2 

192.  DO  372  K =1  , 1 ON 

193.  DO  372  I =1  ,NTT 

194.  WRITE  (6.  25  0 TTT(I),THI(I,K),K 

195.  350  FORMAT(1[X,'TTT=',E12.4,5X,'THI(I,K)=',E12.4,5X,'k=',I3/) 

196.  373  CONTINUE 

| 197.  372  CONTINUE 

198.  C IN  ABOVE  WFAT  WAS  DONE  IS  THAT  'TTT(l)  & THl(l)  ' ARE  MADE  TO  START 

1 9C  . C FROM  1=1.  ORItINALLY  I(TIM)  WAS  SET  = 3 TO  AVOID  THE  FIRST  TWO  COARSE 

2 CO . C TIME  STEPS 

2 C 1 . C NOW  NORMALIZE  PND  PLOT  THE  THEORETICAL  INTENSITIES  STORED  IN  'THI' 
2C2.  C 

2C3.  IF  (FLAG  1. EC. 0)  60  TO  400 

2C4.  CALL  PLO  TT  (THI  ,2,1 ) 

2C5.  GO  TO  5CC 

2C6.  400  CALL  PLO TT ( T H 1 , 2 , 2 ) 

j|  2C?-  C 

2 C8  . C THE  CALCULATIONS  PLOTTING  ARE  OVER  . 

2C9.  C 

J 

2 1 C . 500  STOP 

211.  END 


(j 


1. 

2. 

3. 

4. 

5. 

6. 

7. 

8. 
9. 

10. 

11. 

12. 

13. 

14. 

15. 

16. 

17. 

18. 

19. 

20. 
21. 
22. 

23. 

24. 
25  . 
26. 

27. 

28. 
29. 
3C  • 

31. 

32. 

33. 

34. 

35. 


SUBROUTINE  S 0L VE C N S T A R T , N EM D ) 

C 

C INTEGRATION  PROGRAM  FOR  THE  COUPLED  1-ORDER  DIFFERENTIAL  EQUATIONS 
C NUMERICAL  INTEGRATION  DONE  BY  " PREDICTOR  CORRECTOR  “ METHODS  OF  ADA 
C MS  & MOULTON  TYPE.  THIS  PROGRAM  HAS  THE  OPTION  FOR  SELF  STARTING.  IT 
C VARIES  THE  STEP  SIZE-H  TO  OBTAIN  THE  REQUIRED  ACCURACY. 

C 

COMMON  J ,JMAX,ND,H,DTMAXfNDT,KSTEP,KSEND,END,ED,£H,EN,KDOB, 

1 K PAV1 ,KHAV 2,KHAV3, POBl ,D0B2 

COMMON  T (60  ,EPS  1 (60  ) , Y ( 60 , 1 8 ) , Y R ( 6 0 , 1 8 
1 ) ,ANPR IM< 1 8) 

COMMON  S (18) , A LPH A (18)  , D E NS (100) ,TEMP(100) ,TIME(100),SA( 18) , S B ( 1 8 ) 
1 ,SEC18),NAME(8)fKION(18),DNE(60),TEC60),SFN(18),SEX(18> 

COMMON  NTIME,I0N,I0NM1,NA,NAINI,001,YMIN,YYMIN 

COMMON  TH<60,18),TIM,FLAG,F(1fc>,YES(l8>»W<18),M<5,18),HI<5,18) 

! ,FLAG1  ,TTTC6C) ,DIF2,A(18) 

C 

INTEGER  C0E1 , D0B2, MXDOB 
INTEGER  END 
C 

C COEFFICIENTS  FOR  THE  PREDICTOR  CORRECTOR  FORMULAE. 

C 

MXDOB  =2**2  5 
CF1  = 2./3. 

C F 2 = 4 . / 3 . 

CF3  = 1 ./8. 

C F 4 = 1 ,/12. 

C F 5 = 5 . / 1 2 . 

CF6  =23.  / 1 2 . 

CF7  = 1./24. 

C F 8 = 5 . / 24  . 

CF9  = 9.Z24. 

C F 1 0=  19. /24. 

NST AR  T = NSTART 
2 GO  TO  <5  ,6,4) , NSTART 


21  - 


c 

C SET  THE  VALUES  AT  THE  THREE  INITIAL  POINTS. 

C 

A H=T(2  ) — T ( 1 ) 

IFCH  101*100*101 
ICO  N E N 0 = d 

i»  R I TE  (6,  dCC) 

6CG  FORMAT  (1  >,  EXIT  DUE  TO  T 1 1*;  ESTEP  H AT  »4  OR  U 101  ') 
RETURN 

101  1F(ABS(T(3)-T(2)-H)-ABS(H)*. 01)102, 102,100 

102  N1 = N D T 
Z 1 = A B S (H  ) 

N D T = G 

105  IF(Z1-A3S(DTMAA))103,104,104 

103  NDT  = N CT  +1 
Z 1 =Z 1 *2  . 

GO  TO  1C5 

104  DOBl=«XD CB/2**NDT 

NDT=N  1 
J = 1 

CALL  RATES 
J = 2 

CALL  RATES 

J = 3 

CALL  RATES 

E P S I ( 2)  = C. 

KSTEP  = 3 
GO  TO  3 
C 

C PROGRAM  STARTS  HERl  IF  NTART=1. 

C 

5 H =■  DTMAX  * 2.  **(-NDT) 

D0B1  = MX  DOB/ Z*  *NDT 
J - 2 


I 


74  . DO  13  N = NA  ,ND 

75.  Y (1 , K)  = Y (2,N) 

.76,  Y«(1 , N)  = YR(2,N) 

77.  Y (3, N ) = Y <2,N) 

78.  13  Y R ( 3 » K ) = Y R ( c » N ) 

7*3 , C 

SC.  C ITERATION  OF  THE  ADAMS  MOULTON  K = ? CORRECTOR. IT  DOES  5 ITERATIONS  HERE 


, SI. 

C 

IF  THE 

CONVERGENCE  IS  NOT  FAST  ENOUGH  TO  MEET  THE  ACCURACY  CHECK, 

I T CU 

■ £2. 

c 

TS  THF 

TIMESTEP  BY  HALF  AND  TRIES  AGAIN  FOR  FIVE  ITERATIONS.  ThIS 

goes  ail 

S3. 

c 

ON  UNTIL  DOB  1 BECOMES  ZERO  THAT  IS  NDT  TIMES. 

fl 

S4. 

c 

11 

S5  . 

c 

86. 

12 

T ( 1 ) = 1(2)  - H 

87. 

T ( 3 ) = T(2)  ♦ h 

. 

88. 

j 

1 T 1 » C 

, 89. 

IT2  =-1 

90. 

24 

ZE  =0. 

91  . 

Z1  = EN  * H 

92. 

CO  18  N = N A , N D 

93. 

Z 2 = C F 1 * Y R ( 2 , N ) 

J 94  • 

Y ( 1 , N ) = Y ( 2 » N ) - ( C F 5 * Y R ( 1 , N ) ♦ Z2  - C F 4 * YR(3,N)>  * H 

95. 

ZY  = Y ( 2 , N ) ♦ ( C F 5 * YR(3,N)  + Z 2 - C F 4 * YR(1,N))  * H 

96. 

c 

97  . 

c 

THIS 

"IF"  OMITS  THE  ACCURACY  CHECK  FOR  SMALL  Y(J,N>,  (E.G.  ZERO) 

98. 

r 

99. 

I F ( Y ( 1 , N ) - YYMIN)  18,18,220 

.1  CC. 

220 

Z 2 = ABSCZY  - Y(3,N))/(ABS(Y(2,N))  + ( A c S ( Y R ( 1 , N ) ) + AS  S ( Y R 

1 Cl  . 

1 

( 2 , N ) ) + ABS(YR(I,N) ) )*  Z1) 

-1  C 2. 

ZE  = AM  AX  1 ( Z E , Z 2 ) 

1 C3  . 

18 

Y ( 3 , N ) = ZY 

1 C 4 . 

I F 

(112)  19,20,20 

1 C5. 

C 

1 C6  • 

c 

THE  FOLLOWING  'IF'  IS  ACCURACY  CHECK. 

1 C7. 

c 

• 

1 C8. 

n 

C w 

I F (ZE  - FH/ 16. ) 21 ,21  , 22 

1C9.  6C1  FORMAT ( 1 X,  'Z E = ',  F 1 0. 5/ ) 

in.  22  IT1  = IT1  ♦ 1 

111.  wPITE(6,  COD  IL  - 23  - 

J 


r 


t 1 

1 

, ':‘- 

K H A V 2 = 1 

wRIT(6,  604) 

I Y 1 

- r *5 

- w • 

CO  TO  34 

I1’. 

35 

*.  H A v 3 = K HA  V1  + 1 

1 c4. 

24 

1F(J  -Jl-AX)  3o,37,37 

1 c5. 

C 

1:% 

C K £ S £ T PRINT 

f 57. 

C 

1 5«. 

37 

IF(KPRINT)  26,27,27 

1 59. 

27 

CALL  print 

16'. 

26 

KPRIN  T = 1 

161. 

J1  = JM  A X - 2 

162. 

. J = 1 

1 £T  . 

HO  36  J = J 1 , J M A X 

164. 

T ( N J ) = T ( J ) 

165  . 

EPS  I ( NJ  ) = EPS  I C J ) 

1 66. 

DO  43  N = NA,ND 

1 f . 

Y ( N J , N ) = Y ( J , N ) 

1 65  . 

43 

YR (NJ  ,N  ) = YR (J  ,N) 

1 6 0 . 

36 

NJ  = KJ  ♦ 1 

| 17". 

J = 3 

171. 

C 

1 72. 

C EXTRAPOLATION:  PREDICT  AT  THE  NEXT  TIPE  STEP  WITH  ADAMS  K = 2 PREDICTOR 

C 

1 lu  . 

c 

EXTRA  FCLATION 

175. 

c 

1 76. 

36 

T(J  4 1)  = T(J)  ■*  H 

177  . 

WRITE (6, £05) 

1 7*. 

605 

F0RMATC1X,'  NORMALCASE.  EXTRAPOLATION.') 

179. 

IF  (H-DTMAX)47,48,46 

1£C. 

47 

DOB  2 = D0B 2+D0B1 

i ; i . 

46 

DC  55  N = N A , NO 

1.-2. 

55 

Y(J*1,N)  * Y ( J , N ) ♦ (CF5*YR(J-2,N)-CF2*YR(J-1,N)*CF6*YR (JtN))*H 

1 3. 

J = J ♦ 1 

16. 

CALL  RATES 

18  3. 

C 

1 S6  . 

C INTERPOLATION  : CORRECTION  WITH  ADAMS  K = 3 CORRECTOR. 

1 e '. 

C 

1 £8.  ZE  = C.  C 

1 £9 . WRITE (6, 606) 

1 5 C . 606  FORKATdX,'  INTERPOLATION.') 

191.  Z 1 = EN  * H 

192.  DO  61  N = t.A,ND 

193.  ZY  = Y(  J-1  , N)  + ( C F7* YR  ( J-3  , N) -C F8 * YR ( J-2 ,N)*C FI  0* YR ( J-1  ,N  )*C F9* 

194.  1 YR(J,K)  ) * H 

195.  C THIS  'IF'  OF*  ITS  THE  ACCURACY  CHECK  fOR  SMALL  Y ( E .6  . ZERO) 

196.  IF(Y(JtN>  - YYMIN)  61,61,211 

197.  211  12  = AB S < 2 Y-Y ( J ,N) ) / ( ABS (Z Y) ♦ ( AB S ( YR ( J ,N) ) ♦ A B S ( YR( J -1  ,N  ) )♦ 

198.  1 ABS(Yfi(J-2,N))+ABS(YR(J-3,N)))  * Zl) 


ZE'  = AM  AX1  (ZE  ,Z2) 

61  Y ( J , N ) = ZY 

EPSI(J)  = 0.5  * EPS  I ( J -1 ) ♦ ZE 
KSTEP  = KSTEP  ♦ 1 
WPITE(6,C07)KSTEP,END,J,EPSI(J) 

607  FORMATdX,'  KSTEP=  ', I 5 ,' EN 0 ', F 1 0 . 5 , 'J ', I 5 , 'E PS  I '» F 10 . 5/ ) 

IF  (KSTEP  - KSEND)  62,  C>3,  63 

63  NEND  = 2 
WRITE (6, 71 C) 

71 C FORMAT (1 >,///'  NUMBER  CF  STEPS  (KSTEP)  BECAME  6T  THAN  LIMIT 
1 KSEND') 

RETUR  N 

62  IF  (END)  64,56  ,56 

64  NEND  = 1 
RETUR  A 

56  IF  (KSTEP-4)  49,40,40 

49  IF  (EFSI(J)  - EH/8.)  6,6,23 

C 

C CONTINUE  INTEGRATION  ******************** 

C 

40  IF  (H  - DTMAX)  67,  70,  70 

67  IF  ( K COB)  68,  69,  69 

69  K D OB  = KDOB  - 1 

GO  TO  7 L 


C 

C *71-T0-#79  BELOW  ARE  FOR  D0UBLIN6  THE  STEP  SIZE.  #6-T0-*49  IN  ABOVE 


C FOR  STEP  SIZE  KEPT  SAME.  »80-T0-#95  ARE  TO  CUT  THE  STEP  SIZE  BY  HALF. 

- 26  - • H I — - 


2 26. 

C 

227, 

68 

IF  (EFSICJ)  - ED)  71  , 70, 

70 

228, 

71 

Z SD  OB  1 .A  AD  .DOB 2 

229, 

C 

230. 

C 

.AND.  IS  A KANIPULATIVE  OPERATOR. 

Z IS  ZERO 

AS  D0B1  IS  N0T=D0B2<#47) 

231. 

C 

232. 

IF  (Z  ) 6,72, 

6 

2 33. 

72 

NJ  = J 

234. 

IF  (J  - 4)  73,  73, 

74 

235. 

73 

J = J PAX  - 1 

236. 

KPRINT  = -1 

2 37. 

GO  TO  2 £ 

2 28. 

74 

J = J - 2 

2 39. 

28 

T ( J - 1)  = T(  J ) 

240. 

T ( J)  s T ( NJ  ) 

241. 

EPS  I < J-1)  = EPSl(j) 

242. 

EPS  I ( J)  = EPSI (NJ) 

243. 

DO  79  N = N A , N D 

244  . 

Y (J-  1,N)  = Y ( J ,N> 

245. 

Y(J,N  ) = Y ( N J , N ) 

246. 

YR  ( J - 1 , N ) = Y R ( J , N ) 

247. 

79 

YR(J,A)  = YR(NJ,N) 

248. 

DOB  1 =D0B  1*  2 

249. 

X 

« 

• 

rsj 

ll 

X 

2 50. 

K DOB  = 1 

251. 

KSTEP  « K S T EP -2 

252. 

GO  TO  6 

253. 

C 

254. 

c 

HALF  THE  TIRE  STEP. 

255. 

c 

256. 

70 

IF  (EFSl(J)-EH)  6,6,80 

257. 

80 

IF  (KFAV2)  81,82,82 

2 58. 

81 

NEND  * 5 

259. 

WRITEC6, 730 

260. 

730 

FORMAT (1 >,///*  THE  ACCUMALATED 

ERROR  IS 

GT  THAN  EH  EVEN  IN  THE 

261. 

1 

RECALCULATION  OF  TWO  EARLIER 

TIME  STEPS 

! #> 

262. 

RETUR  A 

263. 

82 

IF  ( K FAVl ) 83,83,84  _ 27 

84  KHAV1  = -1 
KHAV3  = -2 

DPB2=D0E  2-D0B1 
D OB  1 = D 06  1/2 
H = 0 .5  * H 

I F ( DOB  1 . AE  .C)  GO  TO  86 

85  SEND  = 4 
WRITE  (6,  74  0 

740  FORMAT ( 1 X,/// ' EXIT  AT  «85  AS  D0B1  BECAME  ZERO. I. E TOO  MANY  ITERAT 
1I0NS.INPLT  H M A X IS  TOO  LARGE.  ') 

RETUR  A 

83  KHAV2  = -1 

D0B2=D0B  2- DOB  1 *3 
KSTEP  = KSTEP-2 
IF  (KSTEP-2)  96,96,30 
96  KSTEP  = 4 

30  J = J - 2 

IF(J-2)  31,31,86 

31  J = J MA  X - 1 
KPRINT  = -1 

86  T(J)  = T f J - 1) 

EPS  I ( J)  = EPS  I ( J - 1 ) 

DO  32  N = N A , N D 

Y (J  , A)  = Y ( J -1 ,N) 

32  Y R ( J , A ) = Y R ( J - 1 , N ) 

T ( J - 1 ) = TC J)  - H 
EPSI( J - 1 > = Eh 

Z = 0 .5  * H 
DO  95  N=AA,ND 

95  Y ( J - 1 , N ) = (Y(J  ,N)+Y  (J-2,N)-(YRC J , N ) - Y R ( J -2 , N ) V* Z >* 0 . 5 
J = J - 1 
CALL  RATES 
J = J ♦ 1 
60  TO  6 


END 


I 


1.  SUBROUT  1 NE  PATES 

•2.  C 

3.  C RATE  COEFFICIENTS  AND  INCREMENTS  IN  POPULATIONS  YR(J,K)  AT  EACH  TIME 
• 4.  C STEP  J ARE  CALCULATED. 

5.  C 

6.  INTEGER  YES, FL AG1 

7.  COMMON  J ,JNAX,ND,H,DTMAX,NDT ,KSTEP,KSEND,END,ED,EH,EN,KDOB, 

5.  1 K FAV1  ,KHAV2 ,KHAV3 , DOB1 ,DOB2 

9.  COMMON  T (6C),EPSI(60),Y(60,18),YR(60,16 

10.  1 > ,A  NPR  I M ( 1 8 ) 


11.  COMMON  S (18) ,ALPHA(18) ,DENS (100) ,TEMP(100) ,TIME (100), SA( 16) »SB (18) 


1 , SE (18) .NAME (8) ,KI0N ( 18) , DNE (60),TE(60),SFN(18),SEX(18) 

COMMON  NTIME,I0N,10NM1,NA,NAINI,001,YMIN,YYMIN 

COMMON  TH  (60,18), TIM, FLAG, F(18)»VES(16),W(18),M(5, 18), HIC5, 18) 
! ,FLA6 1 , TT  T (60) ,DIF2,A ( 18) 

C 

I F ( T ( J ) ♦ DTMAX  - TIME(NTIME))  6,7,7 
7 END  = -1.0 
C 

C INTERPOLATION  OF  DENSITY  AND  TEMPERATURE 
C 

6 DO  1 L = 1 , NT  I ME 

1 F ( T(J> -TIME(L)  ) 3,2,1 

1 CONTINUE 

2 DN  E ( J ) = DENS (L  ) 

TE (J)  = TEMP(L) 

ABL  = ( DENS(L*1 )-DENS (L-1 ) ) / ( (T IME < L ♦ 1 ) -TIME (L-1) ) *D ENS  (L)  > 
GO  TO  4 

3 DT  = TIME (L)-  T1ME(L-1  ) 

RDT  = (T(J)-TIME(L-I))  /DT 
DIFFD  = DENS(L)-  DENS(L-I) 

DNE ( J ) = DENS (L-1)  ♦ DIFFD*RDT 

TE ( J ) = TEMP(L-I)  ♦ ( TEMP (L ) - TE MP ( L - 1 ) ) *R DT 
ABL  = DIFFD/(DT*DNE(J)  ) 


36. 

37. 
3C  • 
3C  • 

43. 
41  . 
4?  . 
4?  . 

44. 

45. 

46. 

47. 

48. 

49. 

50. 
51  . 

52. 

53. 

54. 

55. 

56. 

57. 

58. 

59. 

60. 
61. 
62. 

63. 

64. 

65. 

66. 

67. 

68. 

69. 

70. 

71. 

72. 


C CALCULATION  OF  I 0 N I Z A T I ON& R E COMB  I N A T 1 0 N RATE  COEFFICIENTS. 

C 

X 5 E = O.C 
XNPRIM  = 0.n 
4 DO  25  K = NA,ION 

I F (FLAG  1. EO • C ) G 0 TO  26 
C 

C FLAG1  « 0 . EVALUATE  SUSE  ECIP  IONIZATION  RATE  - BURGESS&SUMMERS  • 

C 

S(K>  = COEF<K,TE(J),M.(1,K),HI<1,IC),M(2,tO,HI(2,JO,M(3,IC),HI<3,lO, 
1M(4,K),H1(4,K>,M15,K),HI(5,K))*AOO 
GOTO  30 

26  AAA  = 40. 0 * TE(J>  / SE  (K  ) 

IF(AAA-1  .0)  2C,2G,21 

20  AAA  = 1 . C 
C 

C FLA61 =0.EV ALUATE  8 USE  S EM  I EM P E R I C A L IONIZATION  RATE  - LOTZ&KUNZE. 

C 

21  S ( K ) = SA(K)*((ALOG(  AAA  )**3)  ♦ 4Q.G)*SGRT(TE(J))*E 

1 XP(-SE(K)/TE(J))/(SE(K)  ♦ SB(K)*TE(J)) 

30  IF(XSE)  9,9,10 
1C  IF(XNPRIP)  12,12,11 
C 

C EVALUATE  THEEBODT  RECOMBINATION  "BALPHA"  -D  U C H S &GR  I E M , PHY  S . F LU  1 DS 
C ( 1 966 ) VOL  5. 

C 

11  BBB  = XS it  <TE( J>* ( XNPRIM  + 1 . C) * l XNPR IM  + 1 .0)  ) 

IF(BBB-2  5.0  19,19,18 

18  BBB  = 25  .0 

19  BALPHA  = 1 ,4E-31*DNE(J)*( ( XN PR  I M / ( F LO AT ( K ) - 1 .0 ) ) * *6  ) * <<XSE/TE<J> 

1 )**2>*  EXP(bBB) 

GO  TO  13 

12  BALPHA  = 0.0 

C 

C NOW  EVALUATE  DIELECTfiONIC  RECOMBINATION  'DALPHA'  -L AN D I N I & FO S S I . 

C 

IF(YES(X  J.EO.DGO  TO  14 


C Y E S = 1 MEANS  ION  IS  H,HE,NE,K-NI  LIKE  ION  ( R E C OMB 1 N1 N 6 ) 

C 


GO  TO  15 

14  DALPHA  = 1 .4CE-1C*(<  FLOAT  <K))**2)*F  <K)*SQRT  (WOO  >*EXP  (-0.9134* 
:W(K)/TE( J ) )/SQRT ( T E ( J ) * *3  ) 

IN  1 4 F ( K ) IS  NO  OF  ELECTRONS  IN  THE  OUTER  SHELL. K IS  10N.EG.CV. 

GO  TO  13 

IF  VES*1  RECCMBINING  ION  IS  LI-F.NA-AR  LIKE  I0N.F(K)IS  GIVEN  BY 
LONDINIfcFOSS  I PAPER  TABLE  2. 

15  DALPHA  = 1 .6CE-10*((FLOAT(K))**2)*F(K)*EXP(-0.9134*(FLOAT(K)*1.3)*W 
:(K)/(3.3*TE(J))>/SQRT(TE(J)**3) 

EVALUATE  RADIATIVE  RECOMBINATION  "RALPHA** 


13  kALPHA  = 5 .2CE-14*(FL0AT(K)-1.0)*SQRT(XSE/TE(J))*(0.429* 
1 C.5  *AL CG (XSE /TE C J ) ) +0. 469*CBRT  (TE  (J  ) /XSE  ) ) 

ALPHA (K)=BALPHA*DALPHA+RALPHA 
GO  TO  5 

9 ALPHA ( K ) =n  .0 
5 XNPRIM=A APRIM(K) 

WRITE(6,17)K,TE(J) , B ALPHA ,D ALPHA .RALPHA, ALPHA (K) , S <K) 
17  FORMAT (1 CX  ,12 ,6E10.4) 

25  XSE  = S E ( K ) 

TEST  FOR  NEGATIVE  DENSITIES. 

DO  50  K = NA  1 1 ON 

IF  ( Y ( J v K ) ) 51,50,50 
51  Y ( J , K ) = C.O 

50  CONTINUE 

INCREMENTS  IN  POPULATIONS  FOR  USE  IN  S OL VE . . -P RE D I C T 0 R 
-CORRECTOR  METHOD  OF  SOLVING. 


31  - 


c 


1 12. 

113.  IF(NA-I)  7r,70,71 

154.  70  YR<J,1)  = -DNEU)*(Y(J,1>*S<1>-Y<J,2>*ALPHA(2))*ABl*Y<J,1) 

115.  N AA  = 2 

116.  60  TO  72 

117.  71  NAA  = NA 

1 18.  ALPHA ( N A )=C.O 

119.  S<MA  - 1 ) = 0.0 

1 20.  72  00  8 K = N AA  f I ONMl 

« 121.  YR  (J  ,K  )=DNE ( J )*(Y (J  ,K-1  )*S (K-1 >+Y C Jf K*1 ) *ALPHA(K*  1) 

122.  1 -Y(j,K>*<S(K>+ALPHA(K)>)  ♦ A6L*Y(J,K) 

( 123.  8 CONTINIE 

124.  K * 1 0 A 

125.  YR(J,K)=0NE(J)*(Y(J,K-1)*S(K-1)-Y(JfK)*ALPHA(K))  ♦ A8L*YCJ,K) 

126.  RETURN 

' 127.  END 

l 


II 

SUM1 


i I 

I I 


1 . c 

. 2.  C THIS  IS  EVALUATION  OF  E.C.I.P  IONIZATION  RATE  COEFFICIENT  USING 

3.  C SEMI-CLASSICAL  THEORY  OF  A. BURGESS. 

.4.  C 

5.  C THE  FOLLOWING  IS  COMPUTER  FORTRAN  TRANSLATION  DONE  BY  *.qi_AHA 

6.  C OF  THE  CALCLLATION  PROCEDURE  DESCRIBED  BY  SUMMERS  IN  THE  APPELTON 

7.  C LABORATORY  REPORT  I.M.  367  ( 1 974)  PAGES  5 ft  6 . 

S.  C INCORPORATION  OF  THIS  PROGRAM  IN  TO  THE  CODE  FOR  RAPIDLY  IONIZING 

9.  C PLASMAS  WAS  DCNE  bY  RAJU  DATLA. 


1C.  C 

11.  FUNCTION  CCEF(K,TT,I  ,CHI1,J  ,CHI2,L  .CHI', 

12.  1 M ,CH  1 4 .NA.CHI5) 

13.  DIMENSION  X( 10 ) ,OM ( 10) .TEMP ( 1C) , N (5 ) ,CHI ( 5 ) 

14.  C 

15.  C "XM&"OM"  ARE  "X  <K  ),,&,,OMEGA  (K  )M  IN  PAGE  5. THEY  ARE  THE  NODES  AND 


16.  C WEIGHTS  OF  AN  N-POINT  G A U S S - LA G UE R R E QUADRATURE. 

17.  C " N ( I ) K CHI  (I  ) " ARE  NO  OF  EQUIVALENT  ELECTRONS  IN  "I"  SHELLS 

18.  C IONIZATION  ENERGY  OF  AN  ELECTRON  IN  THAT  SHELL. 

19.  C 

2C • DATA  X/C. 137793,0.729455, 1.8 08 343, 3. 4C1434, 5. 55 24 96  ,8. 33 0153, 

21.  1 11.847755,16.279258,21.996566,29.920697/,  OM/0.3C8441, 

22.  2 0.401  12C, 0. 218068, 0. 62C875E-01  ,0  .95 01 5 2 E -02 , 0. 75 3008E -03  , 

23.  3 0.282592e-04,0.424931E-06,0.183956E-08,0.991183E-12/ 

* 24.  C 
25 . C 

• 26.  Z-K 

27.  C 

28.  C 

29.  N ( 1 ) = I 

3C . N(2)=J 

31.  N ( 3 ) = L 

32.  N (4 ) =* 

33.  N ( 5 ) =N  A 

34.  CHI(1)*CH1 

35.  CHI <2>*C  FI  2 - ri  - 


Ti 


j 

i 


I 


36. 

37. 

•23. 

39. 

40. 
41  . 

42. 

43. 
44  . 

45. 

46. 

47. 
43  . 
43. 

5C. 

51  . 

52. 

53. 

54. 

55. 

56. 

57. 

58. 

59. 

60. 
61  . 
62. 

63. 

64. 
65  . 
66. 

67. 

68. 

69. 

70. 

71 . 

72. 

73. 


CHI (3)  = C H 2 
CHI  ( 4 > = C H4 
CHI  (5  ) =C  H5 
DO  7 1 = 1 ,5 

IF  ( N ( I ) .E  U.  C ) GO  TO  7 
CHI (I )=C  H (1 ) / 13.606 
7 CONTINUE 
I T = 1 

TEMPC1 ) - TT 

coef=:.: 

DO  11  NE  =1  ,5 

IF  (N(NE  ).EQ.9)  GO  TO  11 
SUM  =0  .0 

AN=Z/SQR  TC  CHI(NE)) 

C 

C "AN"  IS  "NUr  IN  PAGE  6. 

C 

DO  12  K«  1,  10 

EPS  = TEMP(IT)*XU)/13.606*CHI(NE) 
C 

C 2 E P S ” IS  SYMBOL  EPSILON  IN  PAGE  5. 

C 


U=(EPS-C H (NE) )/CHI(NE) 

C1=E-(B  + 1)*AL0fa(B+1)/(3*2) 

C 2 = 1 .0-1  .0/(B+1)**3 

Q = C1 / ( B*  2) ♦0.65343*C2*YP(Z, A N ,B ) /AN 


C 

C "VP 
C AS 
C 


12 


1 1 


" IS  THE  FUNCTION  Y(DELTA)  IN  PAGE  6. THIS  IS  EVALUATED  IN  SUM  2 

A separate  function. 

6=4*0/ (( E*1)*CHI(NE)**2) 

SUM=SU«*EPS*Q*OM(K) 

continue 

A=13.6:6*CHI(NE)/TEMP(IT) 

COE  F = COE  f*N(NE)*SUM*EXF(-A) 

CONTINUE 

COEF=21.72C?E-G9*COEF*SQRT(13.606/TEMP(IT)> 


C 


SUM2 


C 

C 

C 

C 

c 

c 

c 

c 


c 

c 

c 


"YP"  IS  A FUNCTION  REFERENCED  IN  COEF.  SO 
E.C.I.P  IONIZATION  RATE  COEFFICIENT  EVALUATION  CONTINUED. 

FUNCTION  YP(Z,AN,B) 

DIMENSION  AA (2u> ,BB(2D) 

THE  FOLLOWING  IS  DATA  FOR  A A SB  B WHICH  ARE  A(J)  & 8(J)  IN  PAGE  6 
AND  TABLE  IS. 

DATA  AA/2.  39  16, 1.6742, 1.2583, 0.97 38  , 0.7656, 0.60 78,0  . 4 8 56 1,0. 38976, 

1 0.313E8, 0.25342, 0.20501, 0.16610, 0.13476, 0.10944, 0.08896, 

2 0.072  37  ,0.0 58903  ,. 0479 71  , 0 .03  90S  6 , 0 . 03 1 S 60 / , BB / 1 . 0091  , 

4 0.301 5,0.1314,0.0763,0.0504,0.03561,0.02634,0.01997,0.01542, 

5 0.012  C5  ,0.9  5 3E-02,0.7  5 7E-02,0.602E-02,0.484E-02,0.389E-Q2, 

6 0.312  ZE -02, 0.25 3 5 E-02 ,0.2  04  7 E -02 ,0.1 6 59 E -02,0. 13  44 E -0  2/ 

A=SQRT(B  41  .0) 

D=(Z/AN)  *( (5*AN **2+1 )/(4*Z)+2*AN**2*A/(Z*Z*(B+2)))/(A+ SORT  CB)) 

"0"  IS  "DELTA”  IN  PAGE  6. 


IF  <D ,LE  .0.1 ) GO  TO  1 
I F (D .GE  .2.0)  GO  TO  2 
I =1 0*  D 
R=1C*D-1 
S=I+1-10  *D 

YP  = R*AA  < I* 1)*S*AA(I)-M  (R**i-R)*BB (1  + 1 )* (S**Z-S)*BB ( I)  > /6 .0 
RETURN 

1 A = ALOG(1  .1  229/d) 

VP*A+0.2  5*0*0* (1 .0-2*A*A) 

RETURN 

2 rP=1.5707963*EXP(-2*D)*(1.0+0.25/D-0.09375/D**2+0.0703125/D**3 

1 -0.0472/0**4) 

RETURN 


1 . 
2. 
3. 

4 . 

5 . 
d . 

7. 

8. 
<3. 

10. 
11  . 
12. 

13. 

14. 

15. 
Id  . 
17. 
18  . 

19. 

20. 
21  . 
22. 

23. 

24. 
25  . 
26. 

27. 

28. 

29. 

30. 
31  . 

32. 

33. 

34. 

35. 


SUBROUTINE  PRINT 
C 

c subroutine  tc  print  the  theorical  intensities  and  ionic  populations 
C COMPUTED  by  'solve'. 

c 

INTEGER  TIP 

COMMON  J ,jkax,nd,h,dtmax,ndt,kstep,ksend,end,ed,eh,en,kdob, 

1 K PAV1 ,KHAV2,KHAV3,DOB1 ,D0B2 

common  T (60  tEPSI  (60)  f Y (60 , 1 8)  , YR  ( 60, 18 
1 ) ,ANPRIM(18) 

COMMON  S<18),ALPHA(15),DENS(100),TEMP(100),TIME(100)tSA(18>,SB(18) 
1 ,SE(18),NAME(8),KI0N(18),DNE(60),TE(60),SFN(1S),SEX(18> 

COMMON  ntime,ion,ionmi,na,naini,oo1,ymin,yymin 

COMMON  TH(6  0,18),TIM,FLAG,F(16),YES(18),W(18),M(5,18),HI(5,18) 

! , FL AG 1,TTT (60) , DI  F2 , A ( IE) 

DIMENSION  CONC(1s),XINT(1S),ATOTAL(100) 

C 

1 F < ION  - 9)  1,1  ,2 
C 

1 WR I TE ( 6 * 201)  (KION(L) ,L=1  ,9) 

200  F0RMAT(E12.4,9F11.4,F10.4,3X,I2) 

201  FORMAT (1  H ,1 2H  TIME  (SEC),4X,5H  ION  I2,4X,5H  ION  I2,4X,5H  ION  12, 

1 4 X , 5 H ION  12, 4X , 5 H ION  I2,4X,5H  ION  I2,4X,5H  ION  I2,4X,5H  ION  I 

2 2 , 4 X , 5 H ION  I2,4X,3HSUM,5X,2HNA///) 

301  FORMAT ( 1 M ,1 2H  TIME  (SEC),SX,5h  INT  I2,5X,5h  INT  l2,5X,5H  INT  12, 

1 5 X , 5 H INT  1 2 , 5 X , 5 H INT  I2,5x,5h  INT  I2,5X,5H  INT  I2,5X,5H  INT  I 

2 2 , 5 X , 5 H INT  12///) 

401  FORMAT (1  H ,1 2H  TIME  (SEC),4X,5H  ION  I2,4X,5H  ION  I2,4X,5H  ION  12, 

1 4 X , 5h  ICN  I2,4X,5H  ION  12,4X,5H  ION  I2,4X,5H  ION  I2,4X,5H  ION  I 

2 2 , 5 X , 5 H ION  12///) 

DO  3 L = 4 , J 

01  = 001  /DNE (L) 

TOTAL  = C.C 

DO  6 K = 1 ,9 

CONCOO  = Ol*Y(L,K) 


36. 
77. 
•3e. 

39. 

40. 
41  . 

42. 

43. 

44. 

45. 
46  . 

47. 

48. 

I 49’ 

5C  . 

51  . 

52. 


55. 


| 56. 

57. 

58. 

59. 

eo. 
61. 
62. 

63. 

64. 

65. 

66. 
67. 
63. 

69. 

70. 
71  . 

72. 

73. 


6 TOTAL  = TOTAL  ♦ CONC(K) 

3 WRITE  ( 6 » ZOC)  T (L) , (CONC (K >,K=1 , 9) .TOTAL, NA 
NC  = 1 

NCM  = 10  A 
NCN  = 9 
GO  TO  95 C 

2 URITE(6,401>  (KION(K>,K=1  ,9) 

DO  4 L = 4 , J 

01  = 001  /DUE  (L  ) 

TOTAL  = C.C 

DO  7 K = 1 ,9 

CONC (K)  = 01 *Y  (L  ,K> 

7 TOTAL  = TOTAL  + CONC(K) 

ATOTAL(L)  = TOTAL 

4 WRITE(6,20C)  T(L),(C0NC(K),K=1,9) 

WRITE  (6,  201)  (KION  (K)  ,K  = 1 0, 1 8) 

DO  5 L = 4 ,J 
01  = 001  /ONE (L) 

BTOTAL  = ATOTAL(L) 

DO  8 K = 1 C,  18 
CONC(K)  = C1*Y(L,K) 

3 ETOTAL  = BTOTAL  + CONC(K) 

5 WRITEC6.20C)  T (L > , (CONC (K  ) , K=1 C , 18 >, BTOTAL  ,NA 
NC  = 1 
N CM  = 9 
NCN  = 9 

950  DO  603  K sNC»NCK 

1F(SFN(K  ))  603,603,601 
603  CONTINUE 
CO  TO  60  2 

601  WRITE(6, 301)  (KION (L> ,L=NC, NCM) 

L = 4 

302  IF(L.6T.J)G0  TO  602 
DO  303  K=NC,NCM 
IF(SFN(K  ))  304,304,305 

304  XINT(K)  = C.O 
GO  TO  3C3 

305  V * SEX ( X)  /TE (L)  - 38  - 


■ 


74.  I F ( V -3  .0  7 ) 3 0<  *7  06,337 

75.  3C6  PP  = C • 3 25  *E  I ( V ) 

76 . 00  TO  3C  € 

?7.  307  pp  = 0.2-*0.035/V 

•78.  308  XlNT(K)  = 2.0U-c5*y(L,K)*DNE<L)*SFN<K)*EXP(-V)*PP/SQRT(TE(L)) 

79.  303  CONTINUE 

•8C . C 

81.  C If  THE  FINE  TIME  STEP  'T(L)'  is  equal  to  one  of  the  coarse  time 

82.  C STEP  THEN  FEMEM9ER  IT  IN  'THI'  FOR  LATER  PLOTTING  IN  THE  MAIN  PROGRAM. 

83.  C 

84.  I FCTIMEC TIM) .GT. T(L) )G0  TO  98 


85. 

91 

D1F1=ABS  (TIME(TIM)-T(L)) 

86  . 

T 1=T (L) 

87. 

GO  TO  96 

88. 

98 

IF  (TIM. LT.NTIME)GO  TO  99 

89. 

IF  (L  .EQ  .J  )G0  TO  91 

90. 

99 

OIF  2 = ABS  (TIME(TIM)-T(L)) 

91  . 

DO  100  K =1  , ION 

52. 

100 

THI  (TIM, K) =XINT(K) 

93. 

T2  = T (L) 

94. 

t»  R I T E (6,  305)  T (L)  , (XINT (K > 

55. 

102 

L=L  + 1 

96. 

GO  TO  30  2 

97. 

96 

IF  ( 0 I F 1 .GT.0IF2)  GO  TO  90 

58. 

00  92  K = 1 , ION 

99. 

92 

THI(TIM,K)=XINT(K) 

ICO. 

TTT (TIM)  =T  1 

1 Cl  . 

94 

TIM  = TIM*  1 

i C2  . 

1F(T1M.G1.NTIME)G0  TO  102 

1 C3. 

GO  TO  302 

1 C4. 

90 

TTT  (TIM)  -T2 

1 C 5 . 

GO  TO  94 

1 C6. 

309 

FORMAT ( E 12.4 .5X.9E12.4) 

1C7. 

602 

IF  (IOA-NCN)  550,550,350 

108. 

350 

NC  = 10 

1C9. 

NCM  = ION 

1 10. 

NC  N = 18 

111. 

60  TO  950 

I 


= NC ,NCM  ) 


112. 
1 13  . 
1.14  • 
1 15. 
1 16. 
1 17. 
1 18. 
1 19. 
1 20. 
121. 
1 22. 
1 23. 
1 24. 
1 25. 
1 26. 
1 27. 
1 28. 
1 29. 
1 30. 
131. 
1 32. 
1 33. 
134. 


C 

C AS  TT*E  goes  by  the  lower  ions  will  have  no  populations 

C SO  THOSE  IONS  ARE  NOT  PURSUED  ANY  MORE  IN  SOLVE.  ( #5  50TO#fc2  MAKE  THIS  H / 

c 

550  JA  = J - 2 
NAB  = NA 
DO  59  K = 1 1 1 0 N 
KK=K 

IF(Y(JA,tO  - YMIN)  59,60,60 

59  CONTINUE 

60  NA  = KK 

I F ( NAB  - NA)  7b, 82, 78 

78  NAM  = NA  - 1 
IF(NAM)  £2  ,82, Cl 

81  DO  69  JJ  = 1 , J 

DO  79  KK  = 1 , N AM 
YCJ J,KK)  = 1 .0 

V R ( J J , KK  ) = O.C 

79  CONTINUE 

EPSI  ( J J)  =0.0 
69  CONTINUE 

82  RFTURN 
END 
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PLOT 


1 . 
2. 

3. 

4. 

5. 

6. 
7. 
8 . 
9. 

1C. 

11. 

12. 

13. 

14. 

15. 

16. 

17. 

18. 

19. 

20. 
21  . 
22. 

23. 

24. 

25. 

26. 

27. 

28. 

29. 

30. 

31 . 

32. 

33. 

34. 

35. 


SUBROUTINE  PLOTT(INTE,PFLAG,PFLAGl) 

c 

C ' PL 0 T T"  NORMALIZES  DATA  FOR  INTENSITIES  AND  STORES. SO  TELETYPE 
C PLOTTER  ROUTINES  COULD  BE  USED  FOR  PLOTTING  AT  THE  TELETYPE. 

C 

INTEGER  FFLAG.PFLAG1 

C 

C PFLAG  RPFLAG1  ARE  COUNTERS  TO  KNOW  WHAT  IS  PLOTTED 
C 

REAL  MAX  , I N T E 

COMMON  J ,JMAX,ND,H,DTMAX,NDT,KSTEP,KSEND,END,ED,EH,EN,KDOB, 

1 K FA  VI  ,KHAV2  ,KHAV3 , DOB1  ,DOB2 

COMMON  T (60  , EPS  I (60)  , Y ( 60 , 1 8 ) , Y R ( 6 0 , 1 8 
1 ) ,ANPRIM( 18) 

COMMON  S(16)fALPHA(18),DENS(lOQ),TEMP(lOO),TIME(lOO),SA(18),S8(18) 
1 ,$E(18),NAME(8),<ION(18),DNE(60),TE(60),SFN(18),SEX(18) 

COMMON  NTIME,10N,I0NM1,NA,NAINI,001,YMIN,YYMIN 

COMMON  TH(60f18)tTIM,FLAGtF(18),YES(18),W(18)tM(5,18)tHI(5,18) 

! » FLAG1  ,T  TT  ( 6 C)  , D I F 2 » A ( 1 8) 

DIMENSION  INTE (60,18) 

C 

C 

c 

NT=NTIME -2 
DO  13  K = 1 , I ON 

C FIRST  FIND  INENSITY  MAXIMUM  'MAX'  FOR  THE  ION  AND  USE  IT  FOR 
C NORMALIZATION 
M A X = 0 • 0 

DO  1 1 I - 1 » NT 

11  I F ( INTE ( I,K).GT.MAX)MAX=INTE(1,K) 

C IF  THE  MAX  IS  C.O  SKIP  THAT  ION  FOR  PLOTTING. 

IF(ABS(M AX). LT. IE-10)  GO  TO  13 
DO  12  1=1, NT 
INTE ( I ,K  )=  INTE ( I ,K  ) /MAX 
12  IF(ABS(INTE(I,K)).LT.1E-02)INTE(I,K)=0.0 


37.  WRITE(7)ICN,NTt((TTTCl),INTE(I,K)tI=1,NT)fK=1tI0N>,PFLAG. 


33.  1 ( K I ON  (L)  ,SA(L),SE<L>,A(L)»L  = 1,ION),PFLAG1 

39.  7 RETURN 

.40.  END 

I K 

if 

j I 


; 
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V.  ACCURACY  CHECK 


A simple  system  of  two  ionic  stages  is  considered.  The  density 
and  temperature  are  kept  constant.  Recombination  rates  are  neglected. 
The  coupled  equations  are 


" N1  NI1 


(ID 


N1  NI1  ' N2  NI2 


(12) 


The  solutions  are 


-Nl.t 
e 1 


e'NI2t 


-Nl.t 
e 1 


] 


The  neutral  and  singly  ionized  boron  in  a plasma  of  constant 
14/3 

density  (N  = 5x10  /cm  ) and  temperature  (T  = 250  eV)  are  considered. 
The  ionization  rates  for  these  ions  are  taken  to  be 


.7239  x 10  cm^/sec  and  I„  = .2344 


lft-7  3 
10  cm 


/« 


sec. 


A comparison  of  the  results  from  the  analytical  solutions  given 
above  and  the  numerical  solutions  obtained  with  the  code,  is  shown  in 
Table  1.  The  maximum  difference  is  less  than  one  percent. 
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(PTG) 


□ 


Time  (usee) 

ni/n0 

N2/Nq 

Analytical 

Code 

Analytical 

Code 

.025 

.4045 

.4046 

.4986 

.5049 

.05 

.1637 

.1639 

.5735 

.5807 

.075 

.0662 

.0664 

.5095 

.5158 

.1 

.0268 

.0269 

.4131 

.4183 

.125 

.0108 

.0109 

.3216 

.3256 

.15 

.0044 

.0044 

.2453 

.2484 

Table  1.  Comparison  between  analytical  and  numerical 
solutions. 


APPENDIX 


1^  Plotting  on  the  teletype: 

13 

The  program  TVPLOT  together  with  the  subroutines  TTPLT1, 
TTPLTZ  and  TTPLT3  are  listed  below.  The  plotting  is  done  on  the 
teletype  and  only  the  time  histories  of  the  emission  lines  from  the 
1»  t four  ionization  stages  could  be  plotted.  Examples  of  the  plots 
are  given  in  Appendix  II. 
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T VP  LOT 


INTEGER  FFLAG.PFLAGl 
REAL  1NTE.Y1.Y2.Y3.Y4 

DIMENSION  INTE(6G,16),TIMEX(60),Yl(60),Y2(60),Y3(60>,Y4(60), 
*KlON  < 1 8)  ,SA(18),SE(18),A(16) 

READ<7)1CN,NT,(<TIMEX<I),INTE<I,K),I=1,NT>,K=1,I0N>,PFLAG, 
1<KICN<L),SA(L>,SE(L),A(L),L=1,ION),PFLAG1 
DO  6 1 = 1 ,NT 
K=ION 

Y 1 ( I ) = IN  IE  ( I * K ) 

K=K  -1 

Y2(I>  = INTE(I,K> 

K=K-1 

Y 3 ( 1 ) = IN  IE  (I  »«) 

K = K-1 

Y4CDMN1E  U.K) 

CONTINUE 

CALL  TTPLOT<TIMEX,NT.30,5C,Y1,V2,Y3,Y4> 


IF  'PFLAG'  IS  SET  TO  '1  ' THEN  THE  TITLE* 
'EXPERIMENTAL  TIME  EVOLUTION  OF  IONS  IN  A PLASMA- 
SHOULD  BE  WRITTEN 


1 F ( PF  LA6  .GT.  1 ) 60  TO  103 
WRITE  <6,  50 

50  FORMAT (2 X. /IX , 'EXPERIMENT AL  TIME  HISTORIES  OF  EMISSON  LINES  FROM 
ISUCCESSItE  IONIZATION  STAGES  IN  A PLASMA  ') 

GO  TO  50  C 


IF  'PFLAG'  IS  SET  TO 


THEN  THE  TITLE: 


'THEORETICAL  TIME  EVOLUTION  OF  IONS  IN  A PLASMA- 
SHOULD  BE  WRITTEN 


100  IF  (PFLAE  .NE.  2)  GOTO  500 
IF  (PFLAC1.EQ.DG0  TO  300 
WRITE  (6,  200 
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2CG  FORMAT</  12X, 'PREDI CTED  TIME  HISTORIES  OF  EMISSON  LINES  '/  12X, 

1 'LOTZ 8KU K7 E IONIZATION  RATE  COEFFICIENT  WAS  USED  '/15X,4H  ION,5X 
2,2 H A , 1 G * , 7 H E (EV)) 

WRITE(6,21C)(KION(L),SA(L),SE(L),L=1,ION) 

210  FORK  AT  ( 1 fXtl3,lX,El2.4,1X,Fl2.A) 

60  TO  5 0 C 
300  WRITE  (6  , 31 C) 

310  F0RKAT(2  >,//2X, 'PREDICTED  TIKE  HISTORIES  OF  EMISSON  LINES  '/It, 

1 ' E C I P IONIZATION  RATE  COEFFICIENT  ( BUR G E S S 8 S UKK ER S ) WAS  USED 
2'/15X,4M  I ON , 5 X , 2 H a,1Gx,7h  E (EV)) 

WRITE  ( 6 , 320  (KION(L)  , A ( L ) »S  E (L)  , L = 1 » I ON  ) 

320  F0RMAT(13X,13,1X,F1C.4,1X,F12.A) 


TTPLT1 


1 . C 

> 2.  C THIS  SUBROUTINE  SCALES  AND  DRAWS  PRINTER  PLOTS. 

3.  C THE  SCALING  ASSUMES  A TELETYPE  WIDTH  OF  60  LINES  FOR  THE  GRAPH 

4.  C AND  1C  LINES  FOR  THE  AXIS  LAB  EL S --D I F FE R ENT  VALUES  CAN  BE 

5.  C CONJURED  UP  IF  DESIRED. 

6.  C 

7.  C SUEROUTINE  BY  B.K.REID  NOV. 1970 

8.  C 

9.  SUBROUTINE  R P L OT  R ( X , N P T S , LE  N G TH  , W I D TH  , N Y ) 

1C.  DIMENSION  YY ( 1 500) , X (NPTS  ), AX (2) ,LI NE (63) ,M ARK ( 4) ,NMARK( 4) 

11.  DIMENSION  YMIN  (5  ) , YMAX (5 ) »XX  (35  ) , YLB  ( 10) 

12.  INTEGER  VIDTH,WID,H 

13.  WID=MIN( WI DTH.60) 

14.  NAX=WID//1C 

15.  WID  = 10*NAX 

16.  DATA  MARK/ 

17.  DATA  N M A RK  / ' 1 ' , '2  ' , '3  ' , '4  '/ 

18.  DATA  AX/' + '/ 

19.  C 

20.  C ROUND  THE  WIDTH  UP  TO  A MULTIPLE  OF  10  ROWS. 

21  . C 

22.  NAXL= (LE NGTH-1  )//5 

23.  LEN  = 1*5*MXL 
‘24.  C 

25.  C ROUND  THE  ELEMENT  UP  TO  A MULTIPLE  OF  5 COLUMNS 
'26.  C 

27.  (I*1f1,I.GT.5,YMIN(I)=lE37»YMAX(I)=-1E37) 

28.  DO  10  1=1, NPTS 

29.  C 

30.  C WE  SCAN  THRO  LG  H THE  DATA  GETTING  MINIMUM  AND  MAXIMUM  VALUES 

31.  C 

32.  YMAX ( 1 )= PAX(YMAX (1 ) ,X ( 1 ) ) 

33.  YMIN( 1 )= PIN( YM1N (1  ) ,X (I  ) ) 

34.  DO  10  J = 1 , N Y 

35.  YVAL  = Y8Y  (I  ,J ) 
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26. 

37. 

.38. 

39. 

40. 

41. 

42. 

43. 

44. 

45. 

46. 

47. 

48. 

49. 

50. 
51  . 

52. 

53. 

54. 

55. 

56. 

57. 

58. 

59. 

60. 
61. 
62. 

63. 

64. 

65. 

66. 

67. 

68. 

69. 

70. 

71 . 

72. 

73. 


Y*AX(J*1  )=KAX(YMAX(J+1),YVAL> 

Y*I  N ( J *1  ) = MN(YMIN(J  + 1>,YVAL> 

10  CONTINUE 
NS=NY  + 1 
00  101  L -2  , N S 

CALL  TSCALE(YM1N(L),YMAX(L),YLE>,NAX) 

PRINT  895,^ARK(L-1),CYLB(H),H=1,NAX+1) 

899  FORMAT(1 >,Al ,8X,6G10.4) 

C 

C CONVERT  MAXIMUM  TO  SPAN 

c 

101  CONTINUE 

CALL  TSC  ALE(YMIN(1)fYMAX(1)  tXX,NAXL) 

00  100  L =1  ,NS 
YMAX(L)=YMAX(L)-YMIN(L) 

I F ( YMAX ( L)  .EQ . 0)  YMAX(L)=1 
100  CONTINUE 

V Y ( 1 ) * 0 S INITIAL  LINK  IN  COLUMN  ZERO. 

1 X Y = 1 3 INITIAL  STORE  POINTER 
C 

C THE  DATA  HAS  BEEN  SCALED.  WE  NOW  WIFFLE  THROUGH  IT  A SECOND 
C TIME  AND  SET  UP  THE  PLOT  ARRAY 
C 
C 

DO  1001  IX  A = 1 t NPTS 
XV  = X ( IXA  ) 

MY=1+LEN*(XV-YMIN(1))/YMAX(1) 

DO  1001  L=1,NY 

M = WID*(YGY(IXA,L)-YMIN(L+1>)/YMAX(L-H> 

c 

C WE  ARE  GOING  TO  SEARCH  THE  P L OT T E D -P 0 IN T LIST  FOR  A PLACE  TO 
C PUT  THIS  ONE.  THE  DATA  FORMAT  IS: 

C 

c 

c . . . . 

C . SYMBOL  . COLUMN  . ROW  (X-AXIS)  NO.  . FORWARD  LINK  BY  ROW 

C . . . • 

c 


c 

c 


JX  = 1 

LMMY=C"Y  + <M.LS.12)+(L.LS.1S) 


i 

\ 


74. 

75. 
7-6. 
77. 
.78. 
79. 
.80. 
81  . 
£2. 

83. 

84. 
65. 
£6. 
87. 
88  . 
£9. 

90. 

91. 

92. 

93  . 

94  . 
95. 
96  . 
97. 
98  . 
99. 

K0. 
1 Cl  . 
1*C2  . 
1 C 3 . 
1 C4  . 
1 C5 . 
1 C6. 
1 C7. 
1 C8  . 
1 C9 . 
1 1C. 
111. 


1 C02  CONTINUE 

KX=FLD(2 4,  1?,YY(JX)>  a FORWARD  LINK 
IFCKX.EC  .0)  60  TO  1003 

I F(LMMY.  EQ .FLD (Q,24,YY  (KX  ))  ) GO  TO  10C1 
I F( FLDC1 2, 12 ,YY(KX> ) .GT.MY)  60  TO  1003 
J X = K X 

60  TO  10C2 
1003  CONTINUE 
C 

C WE  HAVE  FOUND  THE  CORRECT  SPOT  IN  THE  CHAIN  FOR  THIS 
C ITEM--OPEN  THE  LINK  AND  INSERT  IT. 

C 

I XY=I X Y+ 1 

FLD (0 , 24  ,Y  Y< IX  Y) )=LMMY 
F LD  < 2 4 , 1 2,  YY  ( I XY  ) ) =KX 
FLD(24,1 2 , YY  ( J X ) ) = IX Y 
1001  CONTINUE 
WID  = UID+  1 

CALL  M IM  TPL( YY , I X Y , XX ,LEN  ) 

RETURN 

SU9R0UTI  HE  M I M T°L  ( Y , N X Y , X ,N  X ) 

C 

C THIS  SUBROUTINE  DRAWS  A PLOT  ON  A TELETYPE  FOR  MIMIC 
C 

C Y is  the  array  of  values,  integer,  RANGE  1-WID 
C my  IS  THE  NIML'ER  of  y values 
C X IS  the  ARRAY  OF  X VALUES 

C NX  IS  THE  NUMBER  OF  X VALUES,  I.E.  THE  GRAPH  LENGTH. 

c 

c 

DIMENSION  X ( N X ) 

INTEGER  Y 
RF  AL  X 

DIMENSION  Y ( 1 5 00  ) 
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- 


12  . 

PARAMETER  NMOD=5 

1?. 

P^INT  81 Cn,(AX,H=1,NAX) 

14. 

5100 

FORMAT (1 6X,lH+,6(A6,A4)) 

15  . 

IX=FLD(24,12,Y(1>)  a-  LINK  TO 

16. 

DO  10  1=1,  NX 

17. 

IXA=I//5 

1V 

DO  15  H=  1 * » I D 

19. 

C 

20. 

C CLEAR  WORKING  PUFFER  TO  ZERO  FOP 

21 . 

C 

22. 

15 

LINE(H)= C 

23. 

6C9 

CONTINUE 

24. 

KK=  FLD  (6  ,6  ,Y  ( IX)  ) ♦ 1 

25. 

JJ=FLD(C  ,6  , Y ( I X)  ) 

26. 

MY=  FLD ( 1 2,  12  , Y (I  X ) ) 

27. 

IF(MY.NE.I)  GO  TO  20 

28. 

IF(LINE(  KK  )>  53,  ,60 

29. 

C 

30. 

C I F 

THIS  IS  T FE  FIRST  POINT  HERE, 

31  . 

C 

32. 

LINE(KK)  =J  J 

33. 

GO  TO  21 

34  . 

60 

CONTINUE 

35. 

L INE (KK  ) =-l 

36. 

C 

37,. 

C I F 

THERE  IS  ANOTHER  POINT  HERE, 

38. 

C 

39. 

GO  TO  21 

40. 

5G 

CONTINUE 

41  . 

C 

42. 

C AT 

LEAST  2 CTHER  POINTS  WERE  ON 

43. 

C 

44. 

LINE (KK) =L INE (KK )-1 

45. 

21 

CONTINUE 

46. 

I X=FLD (2  4,  12  » V ( I X ) ) 

47. 

GO  TO  805 

48. 

2: 

CONTINUE 

49. 

DC  40  J= 1, WID  52 

1 


SMALLEST  FIRST 


EACH  L INE 


JUST  STOPE  IT 


RECORD  THAT 


THAT  SPOT 


48 

49 


L=L INE  (J  ) 

M = AE-S  <L  ) 

IF(l.GT.C)  LINE(J)=MARK(L) 

IF(L.LT.C)  L1NE( J)=NMAPK (M) 

IFCL.EQ. C)  L1NE(J)='  ' 

CONTINUE 

IF(MOD(I  -*4  ,NMOD)  .EU.O)  PRINT  P 9 C 0 , X ( I X A > , ( L I N E ( H ) , H = 1 , W I D ) 
IF(rOD(I-»4,NMOD).NE.0)  PRINT  3931,(LINE(H)fH=1,WlD) 
CONTINUE 

PPINT  81  C0f(AXfH  = 1 , N A X ) 

RETURN 

FORMAT (G  14  .4 ,2h  I , 6 1 A 1 ) 

FORMAT (1  IX  1 1 H I ,61A1  ) 


TTPLT2 


1 . 
2. 

3. 

4. 

b . 
6. 

7. 

8. 
9. 

10. 

11. 

12. 

13. 

14. 
15  . 
16. 

17. 

18. 

19. 

20. 
21 . 
22. 

23. 

24. 

25. 

26. 

27. 

28. 

29. 

30. 

31 . 

32. 

33. 

34. 
35  . 


DRIVER  FOR  TELETYPE  PLOT  SUBROUTINE. 

THIS  SECTION  TAKES  THE  CALL,  COUNTS  THE  ORDINATES,  CALLS 
THE  PLOTTING  FUNCTION,  AND  RECEIVES  ARGUMENT  REQUESTS  FOR 
VALUES  FROM  THE  ORDINATES. 

CALLING  SEQUENCE-- 

CALL  TTP LOT  <X,NPTS, LENGTH, W I DT H , Y 1 , Y2  , . . . , YM  ) 

X THE  ARRAY  OF  VALUES  OF  THE  INDEPENDENT  VARIABLE 

NPTS  THE  NUMBER  OF  POINTS  OF  THE  FORM  <X,Y(K>) 

(1 .E  . THE  DIMENSION  OF  X ) 

LENGTH  THE  LENGTH  OF  THE  GRAPH  IN  INTEGER  PRINT  LINES 
WIDTH  THE  WIDTH  OF  THE  GRAPH  IN  COLUMNS 
Y 1 , Y 2 ETC.:  THE  ARRAYS  OF  DEPENDENT  VALUES 


THE  POINTS  (X  ,Y1  ) ,(X,Y2)  , ETC.  NEED  NOT  BE  IN  ANY  ORDER,  NOR 
NEED  THEY  BE  EDITED:  THE  PLOT  SUBROUTINE  WILL  SORT  AND  EDIT 
AUTOMATICALLY.  THE  SCALE  FACTORS  WILL  BE  COMPUTED  BY  THE 
PLOT  PROGRAM  IN  ORDER  THAT  THE  LARGEST  POSSIBLE  PERCENTAGE  OF 
THE  PLOT  AREA  IS  USED  BY  THE  PLOT. 


NO  MORE  THAN  4 ORDINATES  MAY  BE  USED  IN  A GIVEN  CALL 


A X R $ 

S < 1 ) LIT. 


36. 

TTPLOT* 

L » U 

AO ,4  , XII 

• 

COUNT  THE  ARGUMENTS 

37. 

L X 1 , l 

AO , 1 

• 

ONE  AT  A TIME 

38. 

TZ  ,H  1 

Or  AO 

• 

IF  IT'S  NOT  AN  ARGUMENT 

39. 

TP,  X n 

0,  AO 

• 

OR  AN  ALTERNATE  RETURN 

•40. 

JMGI 

AO, $-2 

• 

THEN  STOP.  ELSE  GO. 

41. 

AH 

AO , (-1  ,1  ) 

• 

COMPUTE  RETURN  ADDRESS 

•42. 

S ,H2 

AO  , RTN 

• 

TO  GO  HOME  AGAIN 

43. 

AN  , U 

AO, 5, XII 

• 

COMPUTE  U OF  ARGS 

44. 

S 

AO  , N Y 

• 

SAVE  FOR  A RAINY  DAY 

45. 

LXI  , L 

XI 1 ,4, AO 

0 

LOAD  INCREMENT  FOR  WB  ACK 

VALUE 

46. 

S 

XI 1 ,WB-*-1 

• 

SAVE  INDEX  11  FOR  'NEXT' 

CALL 

47. 

DL 

AC, 0, XII 

• 

GET  X AND  NX 

48. 

D S 

AO, CALL 

49. 

DL 

AO, 2, XII 

• 

GET  LENGTH  AND  WIDTH 

50. 

DS 

AO.CALL+2 

• 

51. 

LMJ 

Xl 1 , RPLOT  R 

• 

CALL  FORTRAN  DRIVE  PROG 

52. 

CALL 

RES 

A 

• 

SPACE  FOR  X,NPTS  .LENGTH  , WIDTH 

53. 

♦ 

NY 

0 

NUMBER  OF  Y VALUES 

54. 

* 

33  ,Wfcj 

0 

WALKBACK 

55. 

RTN 

J 

C 

0 

RETURN  WORD 

56. 

UB 

'TTP LOT' 

57. 

RES 

1. 

58. 

NY 

RES 

1 

0 

59. 

• 

60. 

• 

61 . 

• 

62. 

. FUNCTION  TO  RETURN  N'TH  OBSERVATION  OF  K 'TH  Y 

63. 

0 

64. 

0 

65. 

. CALLING 

SEQ LENCE 

: X = YQY ( N »K  ) 

66. 

• 

67. 

. THE  ALTERNATE  RETURN  IS  TAKEN  IF 

THERE 

IS  NO  SUCH  VALUE 

68. 

• 

69. 

• 

70. 

TOY* 

L 

AO  ,*  1 , XI  1 

• 

K 

71  * 

L 

A 3 , W B ♦ 1 

• 

RESTORE  CALLING  INDEX 

72. 

A , U 

A3 ,3 , AO 

• 

FIND  THE  ARGUMENT 

73. 

L ,H2 

A1  , , A3 

0 

A 1 HAS  ARGUMENT  BASE  ADDR 

- 55 


TTPLT3 


E UP  ROUTINE  TSCALE<VMIN,VMAX,VSCA,NAXIS) 

DIMENSIC  N VSCA  (10) 

LOGICAL  NEG 
PARAMETER  NS  C F = 1 2 
C 
C 

C SUBROUTINE  TC  SCALE  A PLOT. 

C 

C THE  INPUT  VARIABLES 

C V I N AND  V*A  * CONTAIN  THE  OBSERVED  VALUES  OF  M I N 1 »» UM  / M A XI  M UM 
C VALUES  OF  THE  VARIABLE  TO  BE  SCALED.  ON  OUTPUT,  THEY  WILL 
C CONTAIN  THE  ADJUSTED  VALUES  TO  EE  USED  FOR  'NICE'  NUMBERS 
C ALONG  THE  AVIS  OF  THE  PLOT. 

C 

C VSCA  IS  THE  ARRAY  «h I Ch  WILL  CONTAIN  AXIS  LABELS.  NAXIS  IS 
C THE  NUMBER  OF  AXIS  BLOCKS,  SO  THAT  THERE  WILL  BE  'NAXIS+1' 

C ENTRIES  MADE  IN  VSCA. 

C 

DIMENSION  XSCF(NSCF)/0.,1.,1.25,1.5,2.,3.f4.,5.,6.,7.5,8.,10./ 
NL00PO  £ LOOP  PREVENTER 
D O'MN-VMI  N 

D 0* A X = VMA  X 

6 5 j SCF  = (VMAX-VMIN)/NAXIS 

D X 1 =0  . 

ZERO  = VMI  A 

IF(SCF.Eli.C)  GO  TO  6 0 S 
XNC-ALOG  10  (S  C F ) 

INC  = XNC 

IFCXNC.LE.C.)  INC= INC-1 
IF(INC.GT.C)  XNC10=10.**1NC 
IF(INC.L  F.C)  XNC10  = -(13.**AeS(INC)) 

DX  = PSCAL  E( SCF ,-XNC 1u) 

DO  65  1 I 5 r = 1 , N S C F - 1 

651  I F ( DX  .GT  .X  SC F ( ISC  ) ) D X 1 - X S C F ( I S C ♦ 1 ) 


0. 

1. 

2. 

3. 

4. 

5. 

6. 

7. 

8. 
9. 
0. 


C A PROBABLE  S CALe  FACTOR  HAS  BEEN  COMPUTED,  AND  IT  IS  'DXl' 

C WE  NOW  MUST  ASSIGN  A 'ZERO'  TO  THE  SCALE,  AND  SEE  IF  THE 
C ASSIGNMENT  Of  A ROUNDED  ZERO  IS  GOING  TO  CHANGE  THE  SCALE 
C FACTOR. 

C 

IFCVMIN.  EQ  .1)  GO  TO  675 
APMIN=10  C.  OPSCALE  CVMINt-XNClO) 

c 

C SCALE  THE  ZERO  TO  THE  SAME  ORDER  OF  MAGNITUDE  AS  THE 
C INCREMENT,  BIT  MULTIPLY  BY  IOC  SO  THAT  .75  AND  .25  AND  SO 
C ON  MAY  BE  KEFT  IN  INTEGER  CONVERSION. 

C 

J BM I N = AP S( ABKIN) 

NEG=ABM1 K. lt . C 

LAST3=M0 C( JBMIN, 1033)  « GET  LAST  3 DIGITS 
1F(LAST3.LT.C)LAST3  = LAST3-M000 
IDELT  = 0X  1*  100. 0 a GET  3-DIGIT  INCREMENT 
DO  652  I SC  =0 ,8 
LTEST  = ID  ELTMSC 
I FCLTEST  .GT.LAST3)  GO  TO  65  3 
C 

C BACK  OFF  ON  THE  LAST  THREE  DIGITS  UNTIL  WE  COME  TO  THE 

C FIRST  SPOT  MUCH  IS  A MULTIPLE  OF  THE  SCALE  FACTOR,  INCLUDING  0 

C 

652  CONTINUE 

653  CONTINUE 

IF(.NOT.NEG)  LTEST=LTFST-IDELT 
C 

C USE  THE  LARGEST  PREVIOUS  VALUE  FOR  POSITIVE  SCALING,  BUT 
C USE  THIS  VALUE  FOR  NEGATIVE  SCALING. 

C 

JPMIN=JBMIN-LAST3 

ABM  IN  ~(J EM1N*LTEST)/130. 

ZERO*PSC  ALEC ABMIN,X NCI  0) 

VM1N*S1GK(ZER0,VMIN) 

675  CONTINUE 
DX=DX1 

0X1 -NAX1  S*  DX 1 


DX1 *PSCA  LE (0X1 ,XUC10) 
SPAN= VMI  M DX 1 


74. 

75. 
. 76. 

77. 
.7'. 
79. 
f C. 
■ 1 . 
f>  •>  _ 


C THE  SHENANEGANU  ABOVE  ARE  INSERTED  TO  PREVENT  A FLOATING-PT 
C TRUNCATION  H CASES  WHERE  THE  SCALE  fACTOR  IS  LESS  THAN  ONE 
C AND  NOT  EXACTLY  R E P R t S t N T A B L E IN  FLOATING  POINT.  FOR  EXAMPLE. 
C IF  TrtF  EVF.NTLAL  INTERVAL  IS  .1,  THEN  MULTIPLYING  NAXIS  BY 
C 3.1  -ILL  NOT  GIVE  THE  SAME  ANSWER  AS  MULTIPLYING  NAXIS  BY  1 
C AND  1HEN  DIVIDING  bY  13.  THIS  PRECISION  IS  IMPORTANT. 


1 F ( SPAN  . CE  .VKAX)  GO  TO  649 
NLOOP  =NL  COP* 1 
I f ( NLOOP  .LE. 1 ) GO  TO  650 
PRINT  1 Z L 

FORMAT  ('  ERROR  IN  PLOT  SCALING  ROUTINE:") 


89. 

D12! 

FORMAT  ( 1 X,A6,1x,&12.5,013) 

90. 

D 125 

FORMAT  (1 

X , A6 , 1 X , 1 1 2 ) 

91  . 

D 

PRINT 

1? 

2, 'SPAN", SPAN, SPAN 

9?  . 

D 

PRINT 

12 

2, "vmax',omax,omax 

5 3. 

D 

PPT  NT 

12 

2 , " V M I N' ,OMI N , 0M1  N 

94  . 

D 

PRINT 

12 

2, "7ER0',ZER0,ZfcR0 

55. 

D 

PRINT 

12 

2,  "XNC10',XNCl0,XNCl0 

56. 

D 

PF  I N T 

12 

2, "LAST3',LAST3 

57  . 

D 

PRINT 

12 

5,  ' I D E LT  ' , 1 D E L T 

• 

or 

U' 

D 

PRINT 

1 2 

2, "ltest",ltest 

59. 

D 

P r I N T 

12 

2,  "JBMIN",  JBMIN 

ICO. 

D 

P F I N T 

12 

2,  'DX 1 ", DX  1 , DX  1 

1 C 1 . 

649 

CON  T INUt 

1 C2. 

V"A  X = 

SPAN 

1 C3. 

695 

CONTINUE 

1 C4. 

DX  = PSCAL  E(DX.XNCIO) 

1 C^. 

DO  655  1 ► * 1,NAXIS*1 

10A. 

V S C A ( i H ) 

5VIIIM*(1H-  1 > *DX 

1 C 7. 

655 

C ON  T I NU  r 

1 C 8. 

return 

1 C9. 

C 

1 10. 

c 

111. 

C FUNCTION 

TO 

SCALE  BY  EXPONENTIAL 

- ^ 


c 


FUNCTION  PSCALfc(MANTIS, CHARAC) 

RFAL  MANTIS 

15.  1FCCHARA C.GT.O)  P S C A LE  =M  A NT  I S *C  H A R A C 

16.  IF  CCH  AR  A C.  LT  .0)  PSCALE  = MANTIS/ABS(CHARAC) 

117.  RETURN 

118.  END 

119.  END 


APPENDIX  II 


I 

I 

I 


! 


* 


i 


\ 


| 


; 
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Determination  of  ionization  rate  coefficients: 

One  possible  use  of  this  code  as  mentioned  earlier,  is  to  deter- 
mine ionization  rate  coefficients.  This  involves  changing  the  ioni- 
zation rate  coefficients  in  the  coupled  rate  equations  until  the  pre- 
dicted time  histories  match  with  the  experimentally  observed  time 
histories.  What  is  thus  measured  is  therefore  an  effective  ioniza- 
tion rate,  i.e.  a sum  of  the  ionization  rates  from  the  ground  state, 
including  inner-shell  ionization,  and  all  the  populated  excited  states. 

A detailed  description  of  this  method  was  given  in  Ref.  2.  Another 
useful  result  in  matching  the  time  histories  is  that  the  percentage 
composition  of  the  various  ionization  stages  at  any  instant  of  time  would 
also  be  predicted. 

In  using  the  procedure  for  the  determination  of  ionization  rates 
the  following  points  are  to  be  carefully  considered. 

1.)  As  we  are  theoretically  predicting  the  time  histories  we 
should  carefully  evaluate  how  critically  all  these  predictions  depend 
on  initial  conditions  i.e.,  degree  of  ionization,  electron  density  and 
temperature  at  early  times.  As  we  know,  it  is  very  difficult  to  be 
accurate  about  such  initial  conditions.  So  the  earlier  ionization 
stages  may  not  be  considered  for  measurements.  The  ionization  rate 
coefficients  of  these  ions  could  be  adjusted  until  the  decay  of  an 
emission  line  from  at  least  one  ion  prior  to  the  ion  under  considera- 
tion would  match  with  the  experiment.  This  should  also  make  the  early 
part  of  the  emission  from  the  ion  under  consideration  match  with  the 
experiment.  Thus  the  initial  conditions  could  be  fixed  for  the  ions 
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0 


I 


I 


1 


I 


< 


under  consideration. 

2. )  The  electron  density  N and  the  temperature  T should  be  spa- 
tially non-varying  so  that  ionization  proceeds  uniformly  in  the  plas- 
ma under  observation  and  the  modelling  of  the  plasma  in  Eq.  (3)  would 
be  perfect. 

3. )  Also,  the  electron  density  and  temperature  at  the  time  of 
occurrence  of  the  ions  under  consideration  should  ideally  be  non-vary- 
ing in  time  for  a good  determination  of  ionization  rates.'*'  However, 
some  variations  of  N and  T can  be  accounted  for.  As  was  discussed 

in  Ref.  2,  the  first  three  terms  in  Eq.  (7),  Ref.  2 should  not  dom- 
inate the  time  history. 

A.)  It  is  important  to  select  and  obtain  the  time  histories  of 
emission  lines  connecting  to  the  ground  state  which  are  not  self- 
absorbed. 

As  an  example,  the  determination  of  the  ionization  rate  coeffi- 
cient for  C V (He-like  carbon)  is  considered  from  an  experiment  on 
the  15  KJ  0-pinch.^  The  density  and  temperature  profiles  for  this  plas- 
ma condition  obtained  by  the  laser  scattering  technique  are  given  in 
Fig.  1.  Experimentally  observed  time  histories  of  emission  lines 
from  C IV,  C V,  and  C VI  are  given  in  Fig.  2.  The  predicted  time 
histories  that  are  matched  to  the  experimental  ones  are  given  in 
Fig.  3 and  Fig.  4 using  Lotz  and  Kunze  ionization  rate  coefficients 
and  Burgess  and  Summers  ionization  rate  coefficients,  respectively,  which 

yielded  the  same  results.  The  effective  ionization  rate  coeffici- 

-10  3 

ent  thus  determined  was  0.8  x 10  cm  /sec  at  an  electron  tempera- 
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and  semiclassical  methods  after  including  the  effect  of  the  metastable 
3 

state  (ls2p  ( S))  ionization  agrees  with  the  experiment.  The  mat- 


ching of  time  histories  was  sensitive  enough  to  show  a significant  dif 
ference  if  the  C V ionization  rate  is  varied  by  ±25%. 
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DATA  INPUT  TO  THE  COMPUTER 

CAJi&OiIj.  11  MT0R.I 


50 


6 37 

1 

2 

1 

1 

0.  5 3E- 

08 

3* 

00E 

00 

r\ 

C. 

0.  2QE- 

09 

3* 

COE 

00 

3 

0.7 

09 

3. 

00E 

00 

4 

C.  1 82- 

09 

3. 

COE 

00 

5 

5.  OGE- 

10 

3. 

CCE 

00 

6 

2.  G4E- 

10 

3. 

00E 

00 

6 2 

10000  1 

. 02 

+ 00  1 

• G2- 

03 

0.  10E-06 

2.  502 

01 

0. 

002 

0 0 

0. 00E  00 

2.  502 

01 

0. 

50E 

1 5 

0.'  102- 06 

3.  502 

01 

1. 

GOE 

1 5 

0.  20E-06 

3.  502 

01 

1 . 

502 

15 

0 . 302-06 

3.  5QE 

01 

1. 

5C2 

15 

0.  402- 05 

3.  5 02 

01 

1. 

9 GE 

15 

0*  50E- 05 

0.9  02 

02 

1. 

9CE 

15 

0. 502- 05 

1.  052 

02 

1. 

9 GE 

1 5 

0.'7  0E-Q5 

1.162 

02 

1. 

9GE 

15 

0.8  0E-06 

1 • 232 

02 

2. 

CCE 

15 

0.9  0E-C6 

1.  322 

02 

2. 

5GE 

15 

1 . 00E-06 

1.  402 

02 

3. 

GOE 

15 

1.  10E-06 

1 . 6 22 

02 

o • 

5 OE 

15 

1.  20E-06 

1 . 53  2 

02 

4. 

GOE 

15 

1. 30E-06 

2.  202 

02 

4. 

2GE 

15 

1*40  2-  0 5 

2.  202 

02 

4. 

252 

15 

1. 502-06 

2.  202 

02 

4. 

35_ 

1 5 

1.602-06 

2.  202 

02 

4. 

352 

1 5 

1.7  0E-06 

2.  202 

02 

4. 

352 

1 5 

l.a0E-06 

2.  202 

02 

4. 

35E 

15 

1.9  0E-06 

2.  192 

02 

4. 

35E 

15 

2. 00E- 06 

2.  182 

02 

4. 

30E 

15 

2.  10E-06 

2.  172 

02 

4. 

25E 

1 5 

2.  20E-  06 

2.  152 

02 

4. 

25E 

15 

2.  30 E-  06 

2.  142 

02 

4. 

2GE 

1 5 

2. 40E- 06 

2.  1 32 

02 

4. 

15E 

15 

2.50E-  06 

2.  1 12 

02 

4 • 

1 OE 

15 

2.  60E-06 

2.  062 

02 

4. 

C5E 

15 

2.70E-06 

2.  002 

02 

4. 

COE 

15 

2.’  3 0E-  06 

1.S5E 

02 

3. 

95E 

15 

2.9  0E-06 

1.8  22 

02 

3. 

3 OE 

15 

3.  00E-  06 

1. 60E 

02 

3. 

65E 

1 5 

3.  10E-06 

1 . 462 

02 

3. 

50E 

15 

3.  20E-06 

1.  322 

02 

3. 

35E 

15 

3.'  30E-  06 

1.‘  192 

02 

3. 

05E 

15 

3.  40E-  06 

1. 062 

02 

2. 

7 OE 

15 

3. 50E- 06 

0.96E 

02 

2. 

30E 

15- 

1. 13E  01 
C 1 
01 
01 

02 
02 


2.  44E 

4.  79 E 

5.  452 

3.  9 2w 

4.  9C2 


5.GC2  C 7" 
1 « 0 u 00 
1 • 0 C E 00 

1*  0 'J 

1 • Guu  u J 
1 • DOE  CO- 


DATA 

for  Lotz  6> 
Kunze  Ioniza- 
tion rates. 


1 

1 

1 400 


E+0C  2. 8293  ..+  01  2.  129  3 E+G2 


Experimentally  obtained 
electron  temperature  and 
electron  density  profiles 


(P.T.O) 


For  details,  see-read  statements  In  MAIN-R/W. 


... 


• 4u  j 2+  0 CJ  £•  j 29  j 2+  01  i.»  1 29  3 £■•  0 2 

• 1 4 j L.+  00  1 • 24  376 E+  L«  2 l 3 9 7 !_•>•  ; 

• 4u0  E+0'5  2*  47  3 6415+04 

• ^ c 0 — + 0 1 1 • 6 2 47  0 2+  ^ 2 

• 100  _+ 0 1 2«  39  199L+_3 

• 1 0 0 0 1 1 • 20  9 u 2+  u 0 


DATA  for  Burgess  and  Summers 
(E.C.I.P.)  ionization  rates. 


0 1 
0 ~ 

1 1 

19 


O 

05 

V.  0 

■'  n 

v-  / 

09  5 

O 9 

1 : 

J 

17 

39 

22  5 

15 

07 

50 

_ 

O 2 

3 2~> 

5 

£o 

36 

14 

^3 

55 

145 

19 

39 

19 

1 5 

40 

23 

12 

40 

26  5 

10 

39 

30 

0 

30 

33 

6 

37 

35 

5 

35 

36 

4 

5 3 

3^ 

3 

30 

39  5 

2 

265 

40 

1 

235 

40 

1 

20 

4 0 

1 5 

375 

15 

345 

1 3 

3 1 5 

1 1 

27 

9 

2 3 

i Digitized  data  for  the 
experimental  time  his- 
tories of  emission  lines 
from  C IV,  C V and  C VI. 


INTENSITY  (ARB.  UNITS) 


TIME  - /is 


Fig.  2 EXPERIMENTAL  TIME  HISTORIES  OF  EMISSION 
LINES  FROM  SUCCESSIVE  IONIZATION  STAGES 
IN  A PLASMA 
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